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ABSTRACT 


The  subject  of  this  report  is  the  general  problem  of  signal  processing  for  sensor 
arrays.  Under  certain  reasonable  assumptions,  the  model  for  the  noise  covariance 
matrix  of  the  vector  of  array  outputs  is  an  integral  involving  the  spatial-temporal 
power-spectral-density  function.  This  report  examines  the  application  of  this  co- 
variance  model  to  problems  in  adaptive  bcamforming  and  detection. 

A  constant  false  alarm  rate  detector,  based  on  unconstrained  maximum-like¬ 
lihood  techniques,  is  derived  and  analyzed.  Techniques  such  as  this  do  not  fully 
exploit  the  data  model  and  can  show  an  appreciable  loss  in  performance  compared 
to  optimal  techniques. 

The  space  of  noise  covariance  matrices  possible  from  a  particular  array  is 
characterized,  yielding  representations  for  the  space  and  members  of  the  space  in 
terms  of  finite  numbers  ot  spectral  points.  These  representations  are  used  to  derive 
constrained  maximum-likelihood  estimators  that  jointly  estimate  the  parameters  of 
the  density  function. 

Two  approaches  that  use  the  constrained  covariance  estimates  to  perform 
beamforming  are  described  and  compared.  The  loss  in  signal-to-noise  ratio  and 
the  variance  of  the  estimators  are  shown  to  be  less  for  these  approaches  than  for 
those  that  do  not  use  the  covariance  model. 

Detection  methods  based  on  the  generalized  likelihood  ratio  test  and  a  con¬ 
stant  false  alarm  rate  matched-filter  detector  are  analyzed,  and  simulation  results 
are  presented.  The  analysis  shows  that  the  detectors  will  exhibit  some  constant 
false  alarm  rate  behavior.  The  results  show  a  dramatic  improvement  in  detection 
performance  when  compared  to  unconstrained  approaches. 
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PREFACE 


The  material  presented  in  this  report  is  identical  to  that  in  a  dissertation 
submitted  to  Washington  University  -  Sever  Institute  of  Technology.  December 
1990.  in  partial  fulfillment  of  the  degree  of  Doctor  of  Science. 
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1.  INTRODUCTION 


Sensor  arrays  are  commonly  used  in  fields  such  as  radar,  sonar,  communications,  ultrasound 
imaging,  astronomy,  and  geophysics.  There  are  two  basic  problems  wdiich  are  common  to  these 
fields.  The  first  problem  is  to  estimate  the  intensity  of  a  propagating  wave  from  a  known  direction. 
The  estimation  of  the  intensity  is  a  spatial  filtering  problem  commonly  known  as  beamforming.  The 
second  problem  is  to  determine  if  energy  is  being  radiated  or  scattered  from  given  directions.  This 
detection  problem  is  most  common  in  radar  and  sonar  applications. 

The  interest  in  using  sensor  arrays  is  due  to  the  promise  of  increased  performance  compared 
to  that  of  a  single-sensor  system.  A  sensor  array  can  provide  a  large  physical  aperture  needed  to 
produce  the  spatial  resolution  for  these  problems  without  the  mechanical  disadvantages  of  a  single 
large  antenna.  A  sensor  array  can  also  have  the  desirable  feature  of  adjusting  the  spatial  response 
of  the  sensor  system  to  match  a  changing  interference  environment. 

Optimal  beamforming  and  detection  methods  require  knowledge  of  the  noise  statistics  of 
the  environment.  These  statistics  are  seldom  known  and  must  be  estimated.  Beamforming  and 
detection  performance  are  directly  dependent  upon  how  well  these  estimates  are  made.  Our  interest 
is  in  improving  the  current  beamforming  and  detection  methods  by  utilizing  more  of  the  information 
that  is  available  about  the  array  and  its  response  to  the  environment  to  form  the  estimates  of  the 
noise  statistics. 

1.1  Beamforming 

Beamformers  are  typically  designed  to  satisfy  some  optimality  criteria  such  as  maximizing 
the  output  signal-to-noise  ratio,  minimizing  the  expected  value  of  the  squared  error  between  the 
estimate  and  the  true  signal,  and  minimizing  the  output  variance  while  preserving  a  constant  gain 
in  the  desired  direction.  With  a  Gaussian  signal  model,  these  criteria  result  in  a  linear  filtering 
operation  where  the  outputs  of  the  array  are  weighted  and  summed  ;1;.  For  this  reason  beamforming 
is  traditionally  considered  a  linear  filtering  operation.  The  optimal  linear  beamformers  for  each  of 
the  c'  er..'  require  knowledge  of  the  noise  covariance  matrix. 

1.2  Detection 

Detection  is  required  in  sonar  and  radar  applications.  This  is  a  multiple  hypothesis  testing 
problem  where  there  could  be  a  radiation  source  or  scatterer  in  any  location.  Targets  are  often 
assumed  to  be  scarce,  and  it  is  common  to  condition  on  a  particular  location  and  to  perform  a 
binary  hypothesis  test  to  determine  the  presence  or  absence  of  a  target  at  that  location.  When 
the  noise  statistics  are  known  and  assumed  to  be  Gaussian,  this  results  in  the  comparison  of  the 
output  power  of  an  optimal  linear  beamformer  to  a  fixed  threshold.  When  the  noise  statistics  are 
unknown,  the  parameters  of  the  beamformer  are  unknown,  and  adaptive  approaches  to  detection 
must  be  utilized. 
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Several  quantities  are  used  to  quantify  detector  performance.  The  first  of  these  is  the  proba¬ 
bility  of  false  alarm  (PFA),  which  is  the  probability  that  the  presence  of  a  signal  is  indicated  when 
it  is  not  present.  The  second  performance  measure  is  the  probability  of  detection  (PD).  This  is  the 
probability  that  the  detector  indicates  a  signal  is  present  when  it  is.  For  a  particular  detector  these 
probabilities  are  related:  typically,  lowering  the  false  alarm  probability  also  reduces  the  probability 
of  detection  at  a  particular  signal  level. 

We  are  typically  interested  in  optimizing  detection  performance  based  on  the  Ney  man- Pearson 
criterion  [2,3].  The  Ney  man- Pearson  criterion  is  to  maximize  the  probability  of  detection  while 
restricting  the  probability  of  false  alarm  to  be  less  than  or  equal  to  a  fixed  value.  A  likelihood  ratio 
test,  when  it  exists,  satisfies  this  criterion. 

A  desirable  property  of  detectors  is  independence  of  the  probability  of  false  alarm  with  respect 
to  any  unknown  noise  parameters.  A  detector  which  exhibits  this  property  is  considered  to  be  a 
constant  false  alarm  rate  (CFAR)  detector. 

1.3  Parameter  Estimation 

The  parameters  of  the  noise  environment  must  be  known  in  order  to  perform  optimal  beam- 
forming  and  detection.  In  general  the  statistics  of  the  noise  environment  are  not  known  and  change 
with  time.  If  the  weight  vector  used  to  form  the  linear  filter  is  fixed,  then  there  can  be  a  significant 
loss  in  signal-to-noise  ratio  compared  to  the  optimal  weight  vector  for  the  noise  environment.  The 
probability  of  false  alarm  for  a  detector  is  extremely  sensitive  to  the  noise  covariance  matrix,  and 
a  relatively  small  increase  in  the  noise  level  can  result  in  a  large  increase  in  the  PFA.  For  this 
reason  beamformers  and  detectors  must  be  able  to  adapt  to  the  changing  noise  environment.  The 
coefficients  of  the  beamformer  or  the  parameters  of  the  detector  must  be  estimated  from  the  data 
that  are  available. 

It  should  be  emphasized  that  the  only  reason  w’e  would  want  to  estimate  the  noise  parameters 
is  that  they  affect  the  performance  of  the  beamformers  and  detectors.  The  noise  parameters  are 
nuisance  parameters. 

There  are  two  basic  methods  of  adapting  to  the  changing  noise  environment.  The  first  method 
is  to  initialize  the  beamformer  and  then  estimate  the  weight  vector  as  data  are  received.  This  results 
in  the  wreight  vector  being  recursively  updated  for  each  input  data  point.  The  second  method  of 
dealing  with  the  changing  noise  environment  is  to  utilize  all  of  the  data  available  up  to  that  point 
to  derive  the  test  or  beamformer.  I  will  refer  to  methods  that  do  this  as  block  adaptive  as  they 
operate  on  the  entire  block  of  available  data  at  one  time.  The  emphasis  of  this  report  is  on  block 
adaptive  approaches. 

1.4  Outline  of  This  Report 

The  first  three  chapters  of  this  report  are  devoted  to  introductory  material  concerning  the 
data  model  and  current  methods  of  performing  beamforming  and  detection  with  sensor  arrays.  A 
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detection  method  that  is  based  on  current  techniques  is  the  subject  of  Chapter  4.  In  Chapters  5 
through  7,  constrained  covariances  and  estimation  methods  are  discussed.  The  use  of  constrained 
covariance  estimates  to  perform  beamforming  and  detection  and  a  comparison  of  the  performance 
for  these  methods  are  discussed  in  Chapters  8  and  9. 

The  data  model  is  introduced  in  Chapter  2.  This  model  is  based  on  the  interference  environ¬ 
ment  and  the  processing  which  occurs  in  a  matched-filter  receiver.  The  temporal  samples  of  the 
receiver  outputs  to  a  deterministic  signal  are  shown  to  be  array  “direction  vectors.”  A  stochastic 
model  for  the  interference  is  introduced  which  results  in  an  integral  expression  for  the  covariance 
matrix,  involving  the  spatial-temporal  power-spectral  density.  This  expression  can  be  used  to  show 
that  for  some  arrays  the  covariance  matrix  will  exhibit  structure,  and  an  example  is  given.  The 
goals  of  the  beamformers  and  the  detectors  that  are  proposed  in  this  report  are  discussed. 

Current  methods  of  beamforming  and  detection  are  discussed  in  Chapter  3.  The  bearnform- 
ers  discussed  here  have  been  derived  with  the  intent  of  optimizing  some  aspect  of  the  beamformer 
performance.  Adaptive  beamformers  are  discussed.  The  adaptive  detection  methods  that  are  dis¬ 
cussed  in  this  chapter  include  the  “cell-averaging”  approaches  and  approaches  based  on  generalized 
likelihood  ratio  techniques. 

An  adaptive  detection  method  based  on  the  known  covariance  colored-noise  matched-filter 
detector  is  described  in  Chapter  4.  This  method  utilizes  the  unconstrained  maximum-likelihood 
estimate  of  the  noise  covariance  matrix  and  provides  a  CFAR  test.  The  test  statistic  for  this 
approach  is  simpler  than  the  generalized  likelihood  ratio  test  statistic  and  exhibits  similar  detection 
performance. 

In  Chapter  5  the  covariance  structure  introduced  by  the  integral  expression  for  the  covariance 
matrix  is  investigated.  This  provides  methods  of  representing  the  space  of  possible  covariance 
matrices  and  individual  covariance  matrices  that  are  members  of  that  space.  These  representations 
can  be  used  in  practical  implementations  of  covariance  estimators. 

An  estimator  that  restricts  the  possible  space  of  covariance  matrices  to  those  possible  with  a 
given  array  is  developed  in  Chapter  6.  This  estimator  jointly  estimates  the  mean  and  the  covariance 
of  the  Gaussian  density  function.  The  resulting  iterative  estimator  has  the  property  that  the 
likelihood  for  each  iteration  is  a  non-decreasing  sequence,  and  estimates  that  satisfy  the  necessary 
conditions  for  the  maximizer  of  the  likelihood  are  stable  points. 

In  Chapter  7  a  second  estimator  that  jointly  estimates  the  mean  and  structured  covariance 
matrices  is  proposed.  This  estimator  provides  computational  advantages  compared  to  the  estima¬ 
tor  of  the  previous  chapter.  This  estimator  can  require  less  computation  per  iteration,  and  the 
convergence  of  the  estimates  is  more  rapid  for  some  interference  environments. 

Methods  of  utilizing  the  estimation  procedures  of  the  previous  chapters  to  perform  beam¬ 
forming  and  detection  are  proposed  in  Chapter  8.  Two  beamforming  and  two  detection  methods 
are  proposed  and  discussed.  The  detection  methods  discussed  in  this  chapter  have  some  CFAR 
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properties  and  result  in  tests  where  the  performance  is  much  closer  to  that  which  is  provided  by 
the  known  noise  statistics  test  as  well. 

In  Chapter  9  the  beamforming  and  detection  methods  are  compared  by  use  of  computer 
simulations.  It  is  shown  that  a  dramatic  improvement  in  performance  can  be  obtained  by  using 
the  procedures  proposed  in  this  report.  A  summary  of  th*5  results  is  discussed  in  Chapter  10. 
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2.  DATA  MODEL  AND  PROBLEM  STATEMENT 


In  this  chapter,  the  data  model  will  be  derived,  and  formal  problem  statements  will  be  made. 
The  data  model  is  based  on  the  physical  and  electrical  characteristics  of  the  operating  environ¬ 
ment  and  is  similar  in  both  the  radar  detection  problem  and  the  communications  problem.  The 
data  model  for  the  radar  problem  will  be  derived  in  detail,  and  the  modifications  needed  for  the 
communications  data  model  will  be  indicated. 

2.1  Sensor  Array  in  Noise  Field 

It  is  assumed  that  there  is  an  array  of  sensors  arbitrarily  spatially  located  in  an  environment 
with  propagating  electromagnetic  waves.  The  complex  envelope  of  a  propagating  wave  whose 
properties  are  to  be  examined  is  considered  a  signal,  and  all  other  propagating  waves  are  considered 
noise.  The  sensors  will  respond  to  the  propagating  wavefield  yielding  spatial  samples  of  the  field. 
This  is  shown  pictorially  in  Figure  1.  This  report  is  not  concerned  with  different  propagation  modes 
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Figure  1.  .4  sensor  array  with  incident  signals  and  noise. 


of  the  waves;  therefore,  it  is  assumed  that  the  sensors  respond  equally  to  different  polarizations  of 
the  propagating  waves.  The  response  of  the  receivers  with  a  deterministic  signal  will  be  discussed 
initially,  and  then  the  model  that  results  when  noise  is  included  will  be  examined. 

2.2  Deterministic  Signal 

It  will  be  assumed  that  a  pulsed  signal  is  emitted  by  a  radar  transmitter,  reflected  or  scattered 
by  a  single  point  target,  and  then  sensed  by  a  receiver.  The  propagation  medium  is  assumed  to  be 
homogeneous,  and  no  other  reflectors  or  targets  are  present.  The  signal  transmitted  is 
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(1) 


s(t)  = 


vT^{.4(f)e^}  -T/2<t<T/2 

0  otherwise 


.4(f)  is  the  complex  modulation  of  the  transmitted  signal  and  is  assumed  to  be  normalized  so  that 


f“2  !3?{.4(t)e-,u'°f}|2c/f  =  1 
J-T;  2 


(2) 


.4(f)  is  typically  used  to  create  a  transmitted  signal  with  a  large  time  bandwidth  product.  uj0  is 
the  carrier  or  center  frequency,  and  Et  is  the  transmitted  energy. 

The  complex  envelope  of  the  transmitted  signal  will  experience  the  combination  of  the  fol¬ 
lowing  effects  before  being  received  by  sensor  i. 

1.  Signal  attenuation  Gp  and  propagation  delay  rr(f)  and  r,(e).  The  signal  attenua¬ 
tion  is  due  to  path  loss  and  combines  the  loss  for  propagation  from  the  array  to 
the  scatterer  and  then  again  from  the  scatterer  back  to  the  sensor.  Gp  is  typically 
proportional  to  corresponding  to  a  spherical  wavefront  for  both  propagation  direc¬ 
tions.  The  propagation  delay  is  due  to  the  non- zero  time  that  it  takes  for  the  wave  to 
propagate  from  the  array  to  the  scatterer  and  back.  The  scatterer  is  assumed  to  be 
at  a  range  R(t)  =  Ra  +  R  t  for  range  R0  and  relative  velocity  R  .  The  delay  from  the 
array  reference  to  the  target  and  back  to  the  array  reference  will  be  approximated 
as  rT(t)  =  t0  4-  2 FI  t/c  for  target  relative  velocity  R  much  less  than  the  propagation 
velocity  c  and  where  r0  =  2 R0/c.  There  is  an  additional  differential  delay  term  Tj(e) 
due  to  the  translation  of  the  sensor  away  from  the  array  reference.  This  delay  r,(e) 
is  dependent  upon  the  sensor  location  s,  and  the  direction  of  propagation  e.  For  het¬ 
erogeneous  media,  there  would  typically  be  multipath  and  refraction  effects  which 
would  need  to  be  taken  into  account  here. 

2.  Target  reflectivity  peJ° .  It  is  assumed  that  the  target  scatters  the  signal  and  imparts 
a  phase  shift  and  scale  factor  to  the  complex  envelope  of  the  incident  signal. 

3.  Sensor  response  B,(e).  The  sensor  response  is  a  function  of  the  spatial  origination 
of  a  propagating  wave.  It  will  be  assumed  that  the  target  is  in  the  sensor  far  field 
so  that  the  sensor  response  is  a  function  only  of  the  direction  of  propagation.  The 
sensor  is  assumed  to  have  a  bandwidth  that  is  much  wider  than  that  of  the  receiver 
and  the  transmitted  signal;  the  gain  is  then  independent  of  the  temporal  frequency 
of  the  received  waveform. 

Under  the  model  described  above,  the  signal  received  by  sensor  i  can  be  written 


z,(t)  =  K{y/EtGppj*  Bi(e)A{t  -  r0 


-  2Rt/c  -  r,(e))e^°(‘-T°-2R  Uc-nte))} 


(3) 
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A  received  signal  amplitude  term  bz  can  be  defined  combining  the  effects  of  path  loss  Gp,  target 
reflectivity  pejd .  the  transmitted  energy  \J Et  and  part  of  the  propagation  induced  phase  shift  eJ’1'°T°. 
This  received  signal  amplitude  is 

bi  =  y/EtGpp^e-^  (4) 

and  is  a  complex  parameter. 

It  will  be  assumed  that  all  of  the  receivers  are  identical  matched-filter  receivers  and  that  the 
outputs  of  the  receivers  are  time  samples  of  the  in-phase  (I)  and  quadrature  (Q)  components  of 
the  received  signal.  The  matched  filtering  in  the  receiver  and  the  modulation  A(t)  are  designed 
so  that  the  received  signal  will  be  time  compressed  before  it  is  sampled  with  the  result  that  the 
signal  will  appear  in  only  one  time  sample.1  This  time  compression  is  used  to  increase  the  signal- 
energy-to-interference-energy  ratio  in  one  of  the  time  samples.  A  block  diagram  of  a  receiver 
which  can  accomplish  this  is  shown  in  Figure  2.  The  combined  effect  of  the  radio  frequency  (RF), 
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Figure  2.  Simplified  receiver. 


'For  the  purposes  of  this  report,  the  actual  location  of  the  compression  and  the  implementation 
are  unimportant.  The  matched  filtering  that  implements  the  compression  may  be  located  in  the 
RF,  IF,  video  stages  or  digitally  after  sampling  the  video.  In  practice,  the  actual  location  of  the 
compression  may  be  important  due  to  dynamic  range  and  other  implementation  considerations. 
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intermediate  frequency  (IF),  and  the  detector  sections  is  to  implement  the  matched-filter  receiver. 
The  time  compression  of  the  complex  envelope  due  to  target  motion  will  be  assumed  to  have  a 
negligible  effect  in  the  matched  filtering. 

A  block  diagram  illustrating  the  equivalent  processing  accomplished  for  a  single  receiver  is 
shown  in  Figure  3.  All  operations  that  are  performed  after  the  multiplication  by  the  complex 
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where  the  change  in  phase  -w02 R  t/c  has  been  assumed  to  be  negligible  over  the  period  of  the 
modulation.  A(t)  is  assumed  to  have  been  chosen  so  that  Zi[nTs\  —  0  for  nTs  t0. 

If  transmitted  pulses  are  spaced  Tp  seconds  apart,  then  the  received  data  for  range  sample  n 
of  pulse  m  are 


Zi[nTs)  = 


bx B,(e)e~JiAl°T'l'^e~:>nTa“'°2R  /ce-j^omT,2R  /c  nTg  _  T. 
0  nTs  t  t, 


(7) 


Unless  otherwise  noted,  it  will  be  assumed  that  there  is  a  single  pulse:  m  —  0.  A  redefinition  of  6, 
can  now  be  made  to  incorporate  the  phase  shift  ~nTs^02R  jc. 


bi  =  vEtGppeJ'}e-]'J°T°e-]nT’~°2Rl  c 


(8) 


When  there  are  multiple  samples,  the  temporal  samples  corresponding  to  the  same  delay  or 
range  can  be  ordered,  and  these  samples  can  be  further  filtered  into  several  Doppler  “bins.”  The 
methods  presented  in  this  report  can  be  applied  either  to  the  outputs  of  the  Doppler  bins  or  to  the 
temporal  samples  directly. 

2.3  Deteiministic  Signal  for  Multiple  Receivers 

With  the  output  of  a  single  receiver  defined,  the  output  of  the  receivers  for  .V  spatially 
separated  sensors  can  now  be  investigated.  It  will  be  assumed  that  the  target  is  located  in  the 
array  far  field  so  that  the  signal  wavefront  can  be  approximated  by  a  plane  wave. 

If  the  magnitude  of  the  complex  envelope  of  the  received  signal  does  not  vary  over  the  array 
aperture,  then  the  signal  is  described  as  narrowband.  The  signal  will  experience  a  time  delay  with 
respect  to  the  array  reference  dependent  upon  the  relative  position  and  orientation  of  the  sensor, 
the  direction  of  arrival  of  the  wave,  and  the  speed  of  propagation.  In  Euclidean  coordinates,  let 
e  be  a  unit  vector  in  the  direction  of  propagation  of  the  wave,  and  let  s,  be  a  vector  from  the 
array  reference  to  sensor  i.  For  the  signal  received  at  sensor  i,  the  resulting  delay  with  respect 
to  the  array  reference  will  be  -(e  •  s ,)/c  where  indicates  the  standard  Euclidean  dot-product, 
and  c  is  the  propagation  velocity.  Figure  4  illustrates  this  geometrically  for  a  plane  containing  e 
and  s,.  Also  shown  in  this  figure  is  the  wavelength,  Ac  =  27t cjw0,  of  the  traveling  wave.  It  will  be 
convenient  to  express  distances  in  wavelengths  later  in  this  report. 

The  data  samples  from  multiple  receivers  at  time  n  can  be  arranged  in  a  vector.  The  vector 
output  of  the  array  due  to  the  received  signal  will  then  be  a  complex  scalar  times  an  array  steering 
or  direction  vector: 


z[nj  =  6d(u;,  e) 


(9) 
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Figure  4 ■  Plane  wave  propagating  through  array. 


The  array  steering  vector  d(o;,e)  is 

d(u;.  e)  =  [d\  d2  ■  ■  ■  d_\ jT 


(10) 


with  elements 


di  =  Bl(e)e-^oT'^)  =  fli(S)eJ's?8-s*  .  (11) 

Unless  necessary  to  prevent  misunderstanding,  the  dependence  of  the  direction  vector  d(u;,e)  on  ui 
and  e  will  be  suppressed. 

2.4  Stochastic  Sources 

In  addition  to  the  deterministic  signal  which  may  be  present  in  the  environment,  there  will  also 
be  additive  noise  sources  present.  These  additive  noise  sources  can  be  divided  into  two  categories. 
The  first  category  is  noise  that  is  generated  within  the  receiver  due  to  thermal  and  other  effects. 
The  sum  of  the  additive  noise  terms  that  are  generated  within  the  receiver  will  be  denoted  nr[fj. 
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The  second  noise  category  is  noise  that  originates  in  the  environment  and  propagates  as  waves  to 
be  sensed  along  with  any  signals  that  may  be  present.  Natural  background  noise  is  always  present 
at  some  low  level  due  to  cosmic  radiation,  sky  noise,  solar  radiation,  and  blackbody  radiation  of 
the  earth.  There  may  also  be  interference  due  to  other  radio  emissions.  The  contribution  of  the 
spatially  distributed  noise  sources  to  the  receiver  outputs  will  be  denoted  n3[f]. 

The  following  assumptions  will  be  made. 

1.  The  spatially  generated  noise  originates  from  a  zero  mean  stochastic  process  indexed 
by  t  and  e.  Anticipating  the  bandpass  filtering  in  the  receiver,  this  process  is  assumed 
to  be  a  low-pass  process  modulating  a  sinusoid  at  the  carrier  frequency  and  can  be 
written  h{t.  e). 

2.  The  noise  process  is  assumed  to  be  temporally  wide  sense  stationary  and  spatially 
white: 

E{n(ti.ei)n*(f2.e2)}  =  r(f]  -  f2.ei)<5(ei  -  e2)  (12) 

The  spatial-temporal  power-spectral  density  is 

S(a.'.e)=  f  r(r,e)e~J'"'T  dr  .  (13) 

J  X 

3.  It  is  assumed  that  the  sources  of  interference  are  located  in  the  array  far  field.  This 
assumption  is  reasonable  since  systems  designers  typically  attempt  to  minimize  the 
interference  to  a  receiver  system. 

Under  these  assumptions,  the  spatially  generated  noise  that  is  sensed  at  a  receiver  may  be 
written  as  the  integral 


Xi(0  =  5?|^  ^  Bj(e)ii(t  -  Tl(e).e)e]~'olt  r*’€'}dej  (14) 

This  integral  is  over  the  set  of  all  possible  directions  of  arrival  of  the  interference  fie  r,  here  is 
with  respect  to  the  array  reference  and  by  the  far-field  assumption  7-,(e)  =  -£e  •  s,. 

Van  Trees  '4;  has  derived  a  complex  representation  for  bandpass  random  processes  for  a  single 
receiver.  This  representation  and  the  statistics  can  be  utilized  for  multiple  receivers  as  well. 

The  received  noise  will  be  subject  to  the  processing  of  the  receiver.  The  output  of  a  receiver 
z,  due  to  this  noise  is 


=  Zi(t)  nT>  = 


B,{e)h(t  -  r,(e).  e)e  de  *  *h{t) 


I  nT, 


(15) 


and  has  an  expected  value 
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=0 


(16) 


The  cross  covariance  between  the  received  noise  for  two  sensors  can  be  written 
r,  ]  =  E{z,(t);*(f)} 

B,(ei)n(t  -  T,(e1).e1)e~J~°T,iei)  de}  *  *h(t) 


"  E{(1 

Bj(e2)h(t  -  Tj(e2).e2)e^J^oTlie2)  de2  *  *h(t)Sj  | 

which  can  be  shown  to  be 

2 it  JeiQ^  J x 

where  Mldc)  is  the  Fourier  transform  of  the  matched-filter  impulse  response. 

Making  the  substitution  of  variables  a,'  =  (.2  -  -c0)  and  defining 

Slj,\e)  =  4  MU'  -  -c„)  _  -W’-e) 


this  can  be  written 


r,.  -  -  [  r  SU.sme)B:ie\rJ‘-,r'>*'-r‘'*''de<L 

2~  Je-:UB  Jx 


(17) 


(18) 


(19) 


(20) 


The  argument  of  the  exponential  and  the  constant  due  to  the  sensor  directional  response  can 
be  expressed  by  components  of  a  direction  vector  as 


B,(e)BJ(e)fJ~iT>'*,-r'l§!'  =  d.U.ehfJU.e) 


(21 


and  the  co%-ariance  for  the  received  data  vector  can  be  written  as 

R  n  —  E{z  n  z1  n  }  =  —  [  [  SU.  e)dU\  e ) d f  ( .  e)  dcrfe  (22) 

2~  Jecn^  J x 

If  the  assumption  that  the  system  is  narrowband  is  again  made,  the  covariance  matrix  can 
be  written  as 
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(23) 


R  n  =  E{z  n  zf  n}  =  —  f  5(u.'0.  e)d(vC0,  e)dT(^0.  e)  de 

27T  JeeQe 

It  will  be  assumed  that  the  contribution  of  the  spatially  distributed  noise  sources  to  the 
received  data  vector  can  be  modeled  as  a  Gaussian  distributed  random  variable.  This  model  is 
based  on  the  premise  that  the  spatially  distributed  noise  sources  are  Gaussian  distributed  or  that 
there  is  a  sufficient  number  of  independent  and  identically  distributed  noise  sources  that  the  central 
limit  theorem  can  be  applied  to  model  the  resulting  contribution  as  a  Gaussian  source. 

The  receiver-generated  thermal  noise  can  also  be  assumed  to  be  subject  to  the  Gaussian 
distribution  with  mean  0.  The  noise  that  is  generated  in  each  receiver  is  assumed  to  be  independent 
with  a  contribution  to  the  covariance  of  the  sample  vectors  that  will  be  a  diagonal  matrix 

R  ri  =  diagi  erf  .  .  .  <r)v )  ■  (24) 

The  output  noise  which  is  attributed  to  receiver  noise  sources  is  assumed  to  be  independent  of 
the  spatially  distributed  noise  sources.  In  high-performance  receivers  the  noise  statistics  are  very 
tightly  controlled  so  that  these  receiver  noise  variances  may  be  known.  Additionally,  the  sources  of 
noise  in  each  of  the  receivers  may  be  at  the  same  level  with  the  result  that  the  individual  receiver 
outputs  are  identically  distributed.  These  additional  pieces  of  information  can  be  utilized  in  the 
beamforming  and  detection  algorithms. 

ft  is  assumed  that  the  in-phase  and  quadrature  components  of  each  of  the  receiver  outputs 
are  independent  and  identically  distributed.  The  joint  density  function  for  the  real  and  imaginary 
parts  of  the  received  %’ector  will  be  written  using  the  complex  Gaussian  density  function. 

/(z:R.  6d)  =  _ -l_elz_6d),R_1|Z~bdi  (25) 

The  formulation  of  the  complex  density  function  in  this  manner  comes  from  Goodman  5  . 

In  summary,  the  time  samples  at  the  output  of  the  receiver  may  consist  of  the  sum  of  three 
terms: 

z  n  =  b  n  d  -  ry,  n  —  nr  n  (26) 

with  mean 


(  bn  d  nT,  =  rQ 
\  0  nT,  =t  Tn 
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and  with  covariance  given  by  the  following  expression 


Rbij  =  /  S(u;0,e)d(u;0,  e)dt(u;0,e)de  +  diag(cji  . . .  a2v)  (28) 

7een§ 

With  a  suitable  choice  of  the  modulation  .4(f).  the  matched  filtering  in  the  receiver,  and  a 
suitable  choice  of  the  sampling  interval  Ts ,  the  received  data  vectors  can  be  assumed  independent. 

It  is  often  valid  to  assume  that  targets  are  sparse;  therefore,  the  number  of  samples  which  do 
not  contain  a  target  will  be  much  greater  than  the  number  of  samples  which  do.  As  a  simplification 
of  the  data  model,  it  will  be  assumed  that  the  target  signal  energy  wall  appear  in  the  mean  of  only 
one  sample  for  each  transmitted  pulse  and  this  pulse  is  known.  All  other  samples  are  assumed  to 
be  zero  mean.  This  model,  with  the  target  energy  assumed  to  be  indicated  by  the  mean  of  the 
received  data,  is  the  standard  non-fluctuating  model  for  radar  targets  [6,4] . 

A  simplifying  change  of  notation  will  be  made  here,  and  the  temporal  samples  zjn]  will  be 
denoted 


x„  =  z’nj  .  (29) 

The  sampled  data  vectors  will  be  arranged  into  a  data  matrix  Z.  Column  n  of  this  matrix  is  zn. 
The  data  can  be  ordered  so  that  returns  that  (may)  contain  a  target  return.  nTs  =  r0,  are  the 
first  G  vector  samples  corresponding  to  G  transmitted  pulses.  These  data  vectors  will  be  called 
the  primary  data  vectors.  When  there  is  only  a  single  primary  vector,  G  =  1,  and  there  is  little 
chance  of  misunderstanding,  the  primary  vector  will  be  separated  from  the  data  matrix  and  will 
be  denoted  z.  The  next  K  of  the  samples  are  assumed  to  have  zero  mean,  corresponding  to  data 
vectors  that  do  not  contain  a  target  return.  These  data  vectors  will  be  called  the  secondary  data 
vectors. 

The  mean  of  Z  can  be  expressed  as 

B  =  E(Z)  (30) 

The  first  G  columns  of  B  will  be  complex  multiples  of  the  known  direction  vector,  i.e.,  6, d,  and 
the  remaining  K  columns  will  be  zero. 

The  columns  of  Z  are  assumed  to  be  mutually  independent  with  covariance 

Cov(z,)  =  R  (31) 
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2.5  Coordinates 


In  the  above  derivation,  an  explicit  definition  of  the  coordinate  system  has  been  unnecessary. 
For  many  interference  environments,  the  model  lends  itself  well  to  using  spherical  coordinates  for 
the  integral  equation,  as  shown  in  Figure  5.  The  area  differential  de  is  sin  Qd9d<?  and  Equation  (28) 


Figure  5.  Coordinate  system. 


can  be  expressed 

R’n,  =  —  [  S{uj0.  e)d(a;0, 6,  <f>)df(u;0, 9. 6)  sin  9d0do  +  diag(<7?  ...  a\-),  (32) 

where  the  regions  of  integration  have  been  defined  as  ©  and 
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2.6  Communications  Model 


For  a  communications  system,  the  receiver  is  passive,  and  it  is  assumed  that  the  transmitter 
transmits  a  signal  of  the  form 

s(t)  =  y/EtH{b(t)A(t)eiu°t}  (33) 


The  transmitted  information  is  contained  in  the  amplitude  and  phase  of  the  modulation  term  b. 
This  signal  may  change  from  sample  to  sample  according  to  the  transmitted  modulation.  This 
yields  a  model  for  the  received  data  that  is  nearly  the  same  as  that  for  a  radar  system.  In  the 
radar  system,  targets  are  assumed  to  be  sparse,  and  there  will  be  only  a  few  columns  of  Z  that 
are  non-zero  mean.  In  a  communications  system  many  or  all  of  the  columns  of  Z  may  be  non-zero 
mean. 


2.7  Examples 

2.7.1  Spatially  Uniform  Interference 

A  two-sensor  example  will  be  used  in  order  to  illustrate  the  covariance  that  results  from  the 
integral  expression  found  in  Equation  (32)  for  the  covariance  matrix  of  the  spatially  distributed 
noise.  The  sensors  are  assumed  to  have  omnidirectional  gain  which  can  be  arbitrarily  set  to  unity. 
The  interference  is  assumed  to  be  distributed  uniformly  over  all  directions  of  arrival.  A  spherical 
coordinate  system  will  be  utilized,  and  the  coordinate  system  will  be  rotated  such  that  the  sensors 
will  lie  on  the  2  axis,  and  the  term  e  •  (s,  ~  s})  becomes  |s,  -  s_,|  cost?.  This  is  shown  in  Figure  6. 
The  spatial-temporal  spectral-density  function  will  be  denoted  by  S{u>0)  since  it  is  independent  of 
direction.  The  elements  of  Equation  (32)  can  then  be  written 

=  ej^ia,-s,ic9.*sin  OdBdO  .  (34) 


Performing  this  integration  results  in  the  expression 


sin  “  s,  -  s 

ry  =  25(<V - ^ - 


^|s,  -  Sjl 


(35) 


This  then  describes  the  correlations  between  sensors  as  a  function  only  of  the  distance  between 
the  sensors,  the  propagation  time  of  the  medium,  and  u>0.  When  the  distance  between  sensors  is  a 
multiple  of  one-half  wavelength,  |s*  -  s_,|  =  n^,  the  interference  will  be  uncorrelated  between  two 
sensors,  and  the  resulting  covariance  matrix  will  be  a  multiple  of  the  identity  matrix. 
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Figure  6.  Coordinate  system  for  two  sensor  " linear ”  array. 


2.7.2  Toeplitz  Structure 

It  is  commonly  known  that  uniform  linear  sampling  of  a  spatially  homogeneous  noise  field 
results  in  a  Toeplitz  correlation  matrix.  This  is  given  in  Equation  (20),  which  is  repeated  here  for 
convenience. 


(36) 


There  are  two  requirements  that  must  be  satisfied  in  order  to  have  uniform  sampling.  The 
first  is  that  the  antenna  spatial  responses  Bi(e)  must  be  the  same  for  all  i.  Additionally,  with 
suitable  scaling,  the  term  (rj(e)  -  r,(e))  =  e  ■  (s,  —  Sj)  must  be  expressible  as  a  function  of  i  -  j. 
The  resulting  covariance  matrix  has  elements  that  are  a  function  only  of  i  —  j  and  then  has  the 
Toeplitz  structure. 
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The  results  of  this  report  can  be  applied  to  detection  and  beamforming  problems  where  the 
true  covariance  is  assumed  to  be  Toeplitz;  however,  the  intent  of  this  report  is  to  provide  methods 
that  can  be  utilized  for  arbitrary  sensor  geometries  and  spatial  responses. 

2.8  Problem  Statements 

2.8.1  Detection 

The  purpose  of  detection  is  to  determine  if  a  target  is  present  or  not.  Target  presence  is 
indicated  by  a  non-zero  mean  in  a  single  data  sample.  Additional  target-free  data  are  available 
from  a  process  having  the  same  noise  covariance  matrix.  The  direction  toward  a  possible  target  is 
known  and  has  a  corresponding  direction  vector  d. 

This  is  a  binary  hypothesis  testing  problem,  and  the  detector  will  be  designed  with  the  intent 
of  satisfying  the  Neyman- Pearson  criterion.  The  two  hypotheses  and  the  data  that  are  available 
under  each  hypothesis  are  shown  in  Table  1.  This  data  model  will  be  specialized  in  some  of  the 


TABLE  1 


Detection  Data  Model 


Hypothesis 

Data 

Density 

Ho 

No  target  present 

Zl  •  •  -  Z K-G 

N(  0,  R) 

Hi 

Target  present 

Zl-..ZG 

Zc-l  ■  ■  ■  1k~G 

where  R  =  /g€n-  S(e)d(e)d*(e)  de  +  diag[a *  ■  cr^] 

sections  of  this  report.  Restrictions  will  be  placed  on  the  number  of  terms  K  and  G  so  that  a 
detector  based  on  maximum-likelihood  techniques  will  exist. 

It  will  be  assumed  that  the  mean  b  and  covariance  R  are  unknown.  The  true  covariance  is 
assumed  to  be  non-singular.  The  detection  results  for  the  methods  that  are  presented  in  this  report 
will  be  compared  to  the  performance  of  detectors  where  R  is  known. 
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2.8.2  Beamforming 


The  beamforming  problem  is  a  parameter  estimation  problem  where  the  parameter  of  interest 
is  the  signal  mean(s)  6[nJ.  The  direction  of  arrival  of  the  transmitted  waveform  and  the  correspond¬ 
ing  direction  vector  d  are  known.  The  data  that  are  available  are  shown  in  Table  2.  The  covariance 


TABLE  2 

Beamformer  Data  Model 


Data 

Density 

Zj  ...zG 

N(  Mnld,  R) 

ZG-l  •  •  •  Z/C-G 

N(  0,  R) 

where  R  =  Jg_n_  S^d^d^e)  de  +  diag[cr\  ■  a2^-} 

matrix  is  assumed  to  be  unknown  ana  non-singular.  The  estimates  produced  by  the  methods  which 
are  presented  in  this  report  will  be  compared  to  the  estimates  produced  when  R  is  known. 
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3.  CURRENT  METHODS 


3.1  Beamforming 

Often  the  complex  envelope  of  a  propagating  wave  contains  desirable  information  which  is  to 
be  extracted  by  signal  processing.  When  the  environment  consists  of  interference  sources  whose 
temporal  frequency  content  overlaps  that  of  the  desired  signal,  then  temporal  filtering  alone  cannot 
always  be  used  to  separate  the  interference  from  the  desired  signal.  If  the  interference  and  the 
desired  signal  arrive  at  the  sensors  from  different  spatial  angles,  the  angular  separation  can  also  be 
exploited  to  extract  the  desired  signal  from  the  interference.  This  spatial  filtering  process  is  called 
beamforming. 

Spatial  filters  process  data  over  a  non-zero  spatial  aperture.  This  processing  can  be  provided 
by  an  antenna  with  a  continuous  spatial  aperture,  such  as  a  reflector  or  lens  antenna,  or  it  can  be 
provided  by  combining  discrete  spatial  samples  and  beamforming. 

For  many  continuous  aperture  antennas,  the  desired  spatial  response  is  a  pencil  beam  where 
the  response  is  limited  to  a  small  angular  region  about  the  direction  of  the  desired  signal.  This 
response  region  is  known  as  the  main  lobe,  and  other  typically  undesirable  responses  are  known  as 
the  sidelobes.  The  term  “beamforming'’  when  used  with  spatial  arrays  originates  with  this  connec¬ 
tion  to  continuous  aperture  antennas  [7],  This  connotation  of  forming  beams  may  be  misleading 
since  for  “optimum''  beamformers  the  desired  response  may  not  be  a  pencil  beam. 

There  are  many  methods  that  have  been  utilized  for  combining  the  outputs  of  multiple  sen¬ 
sors  to  realize  beamformers.  One  method  is  to  delay  and  sum  the  outputs  of  the  sensors  so  that 
wavefields  propagating  from  a  desired  direction  will  constructively  add.  This  delay  and  sum  beam- 
former  is  applicable  to  either  narrowband  or  wideband  signals.  For  narrowband  wavefields  and 
receivers,  it  may  be  implemented  by  phase  shifting  rather  than  delaying  the  received  signals.  These 
methods  typically  utilize  analog  processing  of  the  sensor  data:  however,  they  are  applicable  to  the 
processing  of  the  receiver  outputs  as  well.  This  report  will  focus  on  methods  of  combining  the 
samples  of  narrowband  receiver  outputs. 

Beamformers  are  typically  a  linear  function  of  the  sensor  array  outputs.  This  equation  can 
be  written  as 


b  =  w+z  ,  (37) 

where  the  weight  vector  w  indicates  the  linear  combination  of  the  array  data  z.  and  b  is  the 
estimate  of  the  complex  envelope.  Often,  conventional  beamformers  for  antenna  arrays  attempt 
to  form  a  pencil  beam  similar  to  continuous  aperture  antennas  and  typically  utilize  a  fixed  weight 
vector.  This  is  used  to  provide  the  resolution  advantages  of  a  large  array  without  incurring  the 
mechanical  disadvantages.  The  performance  of  these  approaches  is  limited  by  the  spatial  response 
of  the  resulting  beamformers.  Interference  that  appears  in  the  sidelobes  of  the  beam  pattern  can 
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mask  the  desired  signal.  Several  approaches,  both  ad  hoc  and  optimal,  have  been  proposed  that 
partially  alleviate  the  sidelobe  problem  at  the  expense  of  main  lobe  width.  Since  these  methods  do 
not  take  the  interference  environment  into  account,  strong  interference  can  still  obscure  the  desired 
signal,  resulting  in  an  unacceptable  error  in  the  estimate.  Performance  criteria  can  be  used  to 
derive  estimators  which  account  for  the  interference. 

One  performance  criterion  which  can  be  used  is  to  form  an  estimate  of  b  that  minimizes  the 
mean-squared  error: 

w0  =  argmin  E{|£>  -  w*z|2},  b  —  w*z  .  (38) 

W 

This  method  was  investigated  by  Wiener  for  scalar  continuous  time  systems  [8].  Levinson  [9] 
reformulated  the  minimum  mean-square  error  problem  for  discrete  systems  and  developed  the 
normal  equations  for  finding  a  transversal  filter  or  colored-noise  matched  filter  that  minimizes  the 
mean-square  error  between  the  estimate  and  the  true  signal.  This  estimator  requires  that  the  noise 
covariance  be  known.  If  the  signal  mean  is  of  the  form  6d  and  the  noise  covariance  is  R,  the 
resulting  optimal  weight  vector  is 

w0  =  A’R_1d  (39) 

with  k  =  1 6 1 2 / ( 1  —  j6j2d*/?-1d)  [7j.  The  resulting  estimate  of  b  is  biased  and  minimizes  the  mean- 
square  error  for  a  given  value  of  b.  The  results  of  Levinson  apply  directly  to  spatial  filtering.  The 
response  of  the  spatial  filter  that  minimizes  the  mean-square  error  will  tend  to  put  nulls  in  the 
magnitude  of  the  response  in  directions  corresponding  to  strong  interference  sources  [7]. 

Another  performance  criterion  is  the  maximization  of  the  signal-to-noise  ratio  at  the  output 
of  a  linear  filter  [10].  Defining  the  expected  values  E\  and  Es+.v  as  the  expected  values  when  only 
noise  is  present  or  when  signal  and  noise  are  present,  the  signal-to-noise  ratio  at  the  output  of  a 
linear  filter  is  [2] 


SNR  = 


|£s*.v{wtz}  -  fvvlw^z}'2 
E.v{|w*z|2} 


(40) 


For  our  problem  where  the  output  is  zero  mean  when  the  input  is  assumed  to  be  noise-only, 
the  signal-to-noise  ratio  is  [10] 


SNR  =  ]6!2 


iwMj2 

w'Rw 


(41) 
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The  linear  filter  which  maximizes  the  signal-to-noise  ratio  is  given  bv  [10] 


w„  =  kR-'d  ,  (42) 

where  k  is  an  arbitrary  constant.  The  Wiener  filter,  which  minimizes  the  mean-square  error,  will 
then  also  maximize  the  output  signal-to-noise  ratio.  For  the  optimal  linear  filter,  the  resulting 
signal-to-noise  ratio  is  [10] 

SNR0  =  |fejadtR_1d  .  (43) 

This  is  the  signal-to-noise  ratio  that  would  be  achievable  if  the  noise  statistics  were  known.  This 
signal-to-noise  ratio  appears  repeatedly  in  detector  performance  analysis,  so  it  is  convenient  to 
define  the  following: 


a  =  SNR0  .  (44) 

Capon  [11]  proposed  a  similar  problem  where  it  is  desired  to  minimize  the  power  out  of  a 
linear  filter  while  constraining  the  gain  to  unity  in  the  desired  reception  direction.  The  resulting 
beamformer  is  known  as  the  minimum- variance  distortionless-response  (MVDR)  beamformer,  and 
the  weight  vector  is  given  by  [11]  2 


R_1d 
W  “  d+R_1d 

The  covariance  matrix  again  must  be  known  to  evaluate  this  expression.  The  variance  of  the 
estimates  produced  by  this  estimator  is  [11] 

E{;wfz  -  E{wfz}|2}  =  (46) 

This  function,  expressed  as  a  function  of  8  (where  d  =  d(0)  )  using  the  estimate  of  the  covariance 
to  be  introduced  in  Equation  (48),  is  known  as  Capons  spectral  estimator  [12]. 


2The  weight  vector  that  is  derived  by  this  method  uses  the  correlation  matrix  R  +  i&]2ddt  in 
this  formula  rather  than  the  covariance  matrix  R.  By  application  of  the  matrix  inversion  lemma, 
(R  +  |6j2ddt)— 1  =  R_1  -  ]b|2R-1ddtR-1/(l  +  |f»|2dtR_ld),  these  two  forms  are  equivalent. 
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The  methods  discussed  to  this  point  have  not  made  use  of  the  density  function  for  the  received 
data.  Using  the  Gaussian  model  derived  in  Chapter  2  and  assuming  that  the  noise  covariance  is 
known,  the  maximum-likelihood  estimate  of  the  mean  term  is  given  by  the  formula  [12] 


dtR_1z 

dfR-M 


(47) 


This  estimate  of  b  also  results  from  using  the  MVDR  beamformer  on  the  data  sample  z.  For  this 
reason  Equation  (46)  is  sometimes  called  the  ML  spectrum  estimator  [12] .  The  Cramer-Rao  bound 
on  the  variance  of  an  unbiased  estimate  of  b  is  given  by  Equation  (46),  and  thus,  the  MVDR 
beamformer  ;s  an  efficient  estimator  of  b. 

Knowledge  of  the  covariance  matrix  has  been  used  in  the  design  of  the  beamformers  discussed 
thus  far.  The  noise  covariance  matrix  is  seldom  known,  leading  to  adaptive  approaches  where  the 
weight  vector  or  noise  covariance  is  estimated  as  well  as  the  mean.  There  have  been  many  methods 
proposed  to  solve  iteratively  for  the  optimal  minimum  mean-square  error  weight  vector  as  the  data 
are  received  [7],  The  performance  of  these  methods  in  terms  of  the  resulting  signal-to- noise  ratio  is 
difficult  to  analyze,  especially  during  adjustment.  For  stochastic  gradient  approaches,  the  rate  of 
convergence  is  dependent  upon  a  step  size  parameter  and  the  distribution  of  eigenvalues  of  the  true 
covariance  matrix.  Additionally,  there  may  be  an  increase  in  the  output  signal-to-noise  ratio  even  as 
the  number  of  samples  applied  for  adaptivity  approaches  infinity.  These  methods  approximate,  in 
an  iterative  fashion,  the  sample  matrix  inversion  method  that  will  be  discussed  next.  The  iterative 
approaches  will  not  be  discussed  further. 

Reed.  Mallett,  and  Brennan  (RMB)[10],  using  the  results  of  Goodman  [5],  proposed  using 
the  maximum-likelihood  estimate  of  the  covariance  matrix  to  form  a  matched-filter  beamformer 
that  adapts  to  the  interference  environment.  This  was  one  of  the  first  papers  to  propose  using 
statistical  methods  based  on  all  of  the  data  that  are  present  to  estimate  the  weight  vector.  This 
method  is  based  on  the  model  where  K  mutually  independent  and  identically  distributed  data 
samples,  assumed  to  be  zero  mean,  are  available  and  which  can  be  utilized  to  form  an  estimate  of 
the  covariance  matrix.  The  covariance  estimate  used  is  the  average  of  the  sum  of  outer  products 
of  noise-only  samples: 


R  = 


(48) 


The  data  samples  are  the  columns  of  Z  in  this  notation.  This  estimate  is  the  unconstrained 
maximum-likelihood  estimate  in  that  it  does  not  incorporate  any  of  the  covariance  modeling  of 
Chapter  2.  This  estimate  is  substituted  into  the  known  covariance  matched  filter  to  form  the 
weight  vector 


w  =  /cR  1  d 


(49) 
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where  k  is  an  arbitrary  constant. 

The  ratio  of  the  signal-to-noise  ratio  at  the  output  of  a  linear  beamformer  to  the  signal-to- 
noise  ratio  of  the  optimal  linear  beamformer  is  given  by  the  equation 


|wfd|2 

c^R^dwtRw 


IdtR-Mj2 

dtR-1ddtR-1RR_1d 


(50) 


When  the  unconstrained  maximum-likelihood  estimate  of  the  covariance  is  used  to  form  the  weight 
vector,  this  signal-to-noise  ratio  loss  factor  is  a  random  variable  whose  statistics  can  be  described 
in  terms  of  the  dimensional  parameters  N  (number  of  sensors)  and  K  (number  of  samples).  The 
loss  factor  in  this  case  is  beta  distributed  over  the  range  [0, 1]  with  a  density  given  by  [10] 


f3(p)  = 


K\ 


(N  -  2)!(A'  +  1  -  A’)! 


(l~p)N-2pK+l-s 


(51) 


with  an  expected  value  of  ■ 

The  original  paper  by  Reed  et  al.  [10]  proposed  using  the  beamformer  to  perform  detection. 
It  has  been  common  to  use  their  suggestion  of  replacing  the  unknown  true  covariance  with  the 
maximum-likelihood  estimate  of  the  covariance  in  adaptive  beamformers  to  estimate  the  complex 
modulation.  Adaptive  beamforming  is  performed  on  non-zero  mean  data  using  a  covariance  matrix 
estimate  based  on  zero  mean  data.  When  this  method  is  used  with  Capon’s  MVDR  beamformer, 
it  can  be  shown  that  the  output  for  a  single  sample  is  the  maximum-likelihood  estimate  of  b  found 
via  the  joint  maximization  of  the  likelihood  function  over  the  mean  and  the  covariance  matrix. 

One  problem  with  this  use  of  the  maximum-likelihood  estimate  of  the  noise  covariance  to 
form  the  weight  vector  is  that  there  can  be  a  large  loss  in  signal-to-noise  ratio  for  a  small  number 
of  vector  samples  forming  the  estimate.  If  the  number  of  vector  samples  is  too  small  (A'  <  Ar), 
then  the  sample  covariance  matrix  will  be  rank  deficient  and  the  inverse  will  not  exist. 

3.2  Detection 

Adaptive  beamforming  and  detection  have  historically  been  combined  due  to  their  close  re¬ 
lationship.  A  common  approach  to  detection  has  been  to  form  a  beam  in  a  particular  direction  to 
get  an  estimate  of  the  mean  in  that  direction  and  then  to  compare  the  magnitude  of  this  estimate 
to  a  fixed  threshold  to  determine  if  the  signal  is  present.  The  probability  of  false  alarm  for  this 
test  will  vary  with  the  noise  covariance,  so  detectors  based  on  this  method  do  not  have  the  con¬ 
stant  false  alarm  rate  (CFAR)  property.  Various  adaptive  approaches  to  setting  the  threshold  have 
been  proposed.  Many  of  these  utilize  some  type  of  “cell-averaging”  or  cell  comparison  where  the 
estimated  signal  power  (|6]2)  for  the  vector  under  test  (the  test  cell)  is  compared  to  a  combination 
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of  the  estimated  signal  powers  for  other  vector  samples  and  for  different  directions  (the  reference 
cells)  [3].  The  resulting  detectors  can  provide  CFAR  behavior.  These  approaches  are  justified  by 
the  observation  that  the  level  of  noise  from  the  reference  cells  would  be  approximately  the  same  as 
that  of  the  cell  under  test.  If  the  beam  that  is  formed  to  estimate  the  mean  does  not  adapt  to  the 
interference  environment,  then  strong  interference  in  the  sidelobes  of  the  array  response  can  result 
in  a  poor  signal-to-noise  ratio  and  poor  detector  performance. 

A  Neyman-Pearson  optimal  detector  is  based  on  the  likelihood  ratio  test  (LRT)  [2,13].  The 
likelihood  ratio  test  is 


A(.)  =  ^  I  7 


Ih0(*) 


(52) 


Ho 


Using  the  discussed  techniques  [2j,  the  test  for  this  problem  can  be  written  as 


Re(b*dtR~1z)  " 
/(IbpdtR-ld)  Ho 


^  7 


(53) 


The  normalization  constants  in  the  denominator  have  been  introduced  to  simplify  the  expressions 
for  the  probability  of  false  alarm  and  the  probability  of  detection.  By  using  various  methods  [10.2], 
the  probability  of  false  alarm  for  this  test  is 

PF A  =  ^erfc  (y)  ,  (54) 

and  the  probability  of  detection  is 

PD  —  ^erfc  (7  -  a)  .  (55) 

a  is  defined  in  Equation  (44),  and  the  complementary  error  function  is  given  by  the  expression 

erfc(i)  =  — =  /  e  1  dt  .  (56) 

v/rr  J  x 

This  detector  is  not  practical  because  the  true  covariance  R  and  the  mean  b  are  rarely  known. 

A  similar  test,  assuming  that  the  mean  is  unknown  and  the  noise  covariance  is  known,  is 
derived  using  the  generalized  likelihood  ratio  test  (GLRT).  The  GLRT  requires  a  maximization  of 
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the  likelihoods  over  the  unknown  parameters  ©  of  the  two  hypotheses  before  forming  the  likelihood 
ratio  ^2,: 


maxe  H\  Ifuiz-Q) 
max©  h0  Ih0(z-  ©)  i 


(57) 


The  test  derived  by  this  method  is 


2 

dtRM 


H  i 

£  *> 


Ho 


(58) 


This  test  is  not  uniformly  most  powerful  because  the  LRT  in  Equation  (53)  can  be  used  assuming 
a  particular  value  for  the  mean,  and  the  probability  of  detection  will  be  higher  when  the  true  mean 
is  near  the  design  value.  The  generalized  likelihood  ratio  test  of  Equation  (58)  can  be  written  as 
the  ratio  of  the  signal  power  estimated  by  the  maximum-likelihood  method  to  the  MVDR  spectral 
estimate 


6  2 

1  d'Rld 


H  i 

^  1 


Ho 


(59) 


and  is  the  test  that  results  when  the  maximum-likelihood  estimate  of  b  is  substituted  without 
the  normalization  in  the  denominator  into  Equation  (53).  This  is  an  example  of  the  relationship 
between  beamforming  (estimation  of  6)  and  detection. 

Other  methods  [10.21  can  be  used  to  determine  the  detection  performance.  The  probability 
of  false  alarm  for  this  test  is 


PFA  =  e-'  .  (60) 

and  the  probability  of  detection  is 

PD  -  Q{\'2a.  v/27)  .  (61) 

where  Q  is  the  Marcum's  Q-function 
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Q(a.b)=  H 
Jb 


ue 


2  2 


(62) 


2  I0(au)du 


The  performance  of  detection  algorithms  that  are  proposed  in  this  report  will  be  compared 
to  the  performance  of  these  known  parameter  detection  algorithms. 

3.2.1  Unknown  Noise  Covariance  GLRT 

Kellv  [14,15’  has  derived  a  generalized  likelihood  ratio  detector  based  on  the  model  that  K 
mutually  independent  zero  mean  secondary  data  vectors  (the  columns  of  Z)  and  a  single,  possibly 
non-zero  mean,  primary  data  vector  are  available.  These  data  are  assumed  to  share  the  same 
covariance  matrix  and  are  subject  to  the  Gaussian  density  function.  The  unconstrained  maximum- 
likelihood  estimate  of  the  covariance  is  used  in  an  adaptive  detection  rule  that  is  given  for  data 
sample  z  and  threshold  7  by  [15] 
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with 


R  =  -pZZt  (64) 

A 

This  detector  requires  an  increase  in  signal-to- noise  ratio  to  achieve  a  given  probability  of 
detection  at  a  fixed  probability  of  false  alarm  when  compared  to  the  known  covariance  detector. 
The  required  signal-tonoise  ratio  increase  is  greater  than  the  loss  due  to  adaptive  beamforming; 
the  additional  loss  is  attributed  to  the  determination  of  the  threshold.  The  loss  is  independent  of 
the  true  covariance,  depending  only  upon  the  dimensional  parameters  and  the  probability  of  false 
alarm.  The  generalized  likelihood  ratio  test  results  in  a  test  in  which  the  probability  of  false  alarm 
is  independent  of  the  level  and  structure  of  the  noise  covariance,  rather  than  just  the  “local"  noise 
level  as  in  the  previous  cell-averaging  techniques. 

We  16,17;  have  derived  an  adaptive  detector  which  is  simplified  compared  to  the  generalized 
likelihood  ratio  detector.  This  decision  rule  is  defined  by  substituting  the  maximum-likelihood 
estimate  of  the  covariance  based  on  the  zc  o  mean  data  samples  in  place  of  the  known  covariance  in 
the  unknown  mean  detector  of  Equation  (58).  The  derivation  and  analysis  of  this  test  are  discussed 
further  in  Chapter  4.  The  detection  rule. 
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compares  the  output  power  of  a  normalized  matched  filter  to  a  threshold.  This  procedure  is 
called  the  adaptive  matched-filter  or  AMF  detector.  It  will  be  shown  that  the  AMF  detector  has 
the  constant  false  alarm  rate  (CFAR)  property  similar  to  the  generalized  likelihood  procedure. 
The  probability  of  detection  of  this  detector  can  show  either  an  increase  or  decrease  in  detector 
performance  compared  to  Kelly's  generalized  likelihood  ratio  detector:  the  difference  in  performance 
depends  upon  the  value  of  the  signal-to-noise  ratio  and  the  dimensional  parameters. 

Both  of  the  detectors  that  use  the  unconstrained  maximum-likelihood  estimate  of  the  co- 
variance  matrix  suffer  from  problems  that  are  similar  to  those  of  the  unconstrained  maximum- 
likelihood-based  beamformers.  A  large  increase  in  signal-to-noise  ratio  may  be  required  in  order  to 
reach  the  same  probability  of  detection  as  the  test  where  the  covariance  is  known.  Furthermore, 
the  tests  will  not  exist  if  K  <  N  since  the  sample  covariance  matrix  is  singular. 

3.3  Covariance  Estimators 

Although  the  covariance  matrix  is  a  nuisance  parameter,  it  plays  an  essential  role  in  the 
detection  and  beamforming  algorithms.  We  will  utilize  the  model  for  the  covariance  derived  in 
Chapter  2  and  estimate  covariance  matrices  that  are  restricted  to  this  model  in  the  hopes  that  by 
estimating  a  matrix  that  is  "closer"  to  the  true  covariance  matrix,  higher  performance  will  result 
compared  to  estimates  that  do  not  make  full  use  of  the  model.  Estimators  based  on  some  of  the 
aspects  of  the  model  in  Chapter  2  have  been  used  previously  for  some  special  arrays  and  interference 
environments:  these  will  be  discussed  in  this  section. 

3.3.1  Unconstrained  Maximum  Likelihood 

The  unconstrained  m  iximum-likelihood  estimate  of  the  covariance  causes  several  problems 
when  used  in  adaptive  beamformers.  The  main  problem  is  that  although  the  expected  loss  in 
signal-to-noise  ratio  may  be  small  for  a  given  number  of  samples  used  to  estimate  the  covariance 
(on  the  order  of  2 N  samples,  the  expected  loss  is  3  dB  [  10] )  the  variance  of  the  loss  is  high.  Any  one 
realization  of  the  adaptive  beamformer  may  have  an  unacceptable  loss  in  signal-to-noise  ratio  [18]. 
The  variance  of  the  estimate  also  generates  an  instability  in  the  spatial  response  of  the  adaptive 
beamformer.  This  instability  has  been  attributed  to  the  sensitivity  of  the  inverse  of  the  estimated 
covariance  matrix  to  the  smallest  eigenvalues  [19.20]  and  is  a  function  of  the  condition  number  of 
the  covariance  matrix. 

The  AMF  and  the  generalized  likelihood  ratio  detectors  also  suffer  from  poor  sidelobe  per¬ 
formance  jl6,21].  If  the  primary  data  vector  has  a  source  of  interference  that  is  not  reflected  in 
the  other  data  samples  either  due  to  non-stationarity  of  the  interference  or  because  the  sample 
contains  an  outlier,  then  detection  performance  will  be  poor.  The  AMF  test  has  poor  rejection  of 
signals  or  interference  appearing  in  the  sidelobes  of  the  range  bin  under  test  but  which  are  not  in 
the  secondary  data  [16] .  The  effect  of  this  is  that  a  strong  outlier  will  cause  the  detector  to  indicate 
that  there  is  a  signal  in  all  directions,  masking  the  detection  of  the  signal  which  might  appear  in 
the  same  vector  sample.  For  Kelly's  GLRT  detector,  a  single  interferer  in  the  primary  vector  will 
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cause  the  test  statistic  to  be  reduced  for  directions  other  than  in  rhe  direction  of  the  strong  signal, 
and  detection  of  the  desired  signal  may  not  be  accomplished  [21].  The  impact  of  the  outliers  is 
reflected  in  the  probabilities  of  detection  and  false  alarm;  non-stationarity  of  the  noise  environment 
is  not  reflected. 

3.3.2  Subspace  Approaches 

An  attempt  to  improve  the  sidelobe  performance  of  the  adaptive  beamformer  over  that  of  the 
unconstrained  maximum-likelihood  approach  is  based  on  modeling  the  interference  as  a  number 
(J  <  N)  of  narrowband  spatial  point  sources.  The  covariance  is  then  given  by  the  rank  J  term 
due  to  the  emitters  plus  a  diagonal  term  attributed  to  receiver  noise.  One  method  that  uses  this 
technique  is  based  on  the  eigenvector  decomposition  of  the  sample  covariance  matrix  [13,22],  This 
method  assumes  that  the  J  largest  eigenvalues  correspond  to  the  emitters  and  that  the  smallest 
N  -  J  eigenvalues  correspond  to  white  noise.  The  smallest  eigenvalues  are  averaged  (or  set  to  a 
known  value)  to  form  the  covariance  estimate,  which  is  then  used  in  an  adaptive  beamformer.  This 
is  the  maximum-likelihood  estimate  of  R  under  this  low-rank  model. 

Another  method  assumes  that  the  emitter  locations  are  known  and  then  estimates  a  covariance 
matrix  which  is  constrained  to  be  of  this  form  [23,24]. 

These  methods  are  sensitive  to  errors  in  the  subspace  chosen  for  the  discrete  emitters  and  in 
the  eigenvalue  spread  of  the  true  covariance. 

3.3.3  Covariance  Modeling 

Covariance  models  can  be  derived  based  on  the  array  geometry  and  the  spatial  structure  of  the 
noise.  A  covariance  estimate  based  on  this  model  can  then  be  utilized  in  an  adaptive  beamformer 
or  detector.  A  common  assumption  is  that  the  interference  is  located  in  the  far  field  and  is  spatially 
independent.  This  assumption  leads  to  a  spatially  stationary  noise  environment.  A  uniform  linear 
array  sampling  this  field  results  in  a  Toeplitz  covariance  matrix.  Estimators  for  covariances  that 
are  constrained  to  have  the  Toeplitz  structure  are  popular  [25.26,27,28].  The  two  predominant 
methods  of  forming  an  estimate  that  has  the  Toeplitz  structure  are  to  use  maximum-likelihood 
methodology  or  to  use  diagonal  averaging. 

Diagonal  averaging  minimizes  the  Frobenious  norm  of  the  difference  between  the  sample 
covariance  and  the  estimate.  An  unbiased  estimate  may  result,  but  it  will  not  always  be  positive- 
definite.  Modifications  have  been  introduced  to  ensure  a  positive-definite  estimate. 

The  Toeplitz  constrained  maximum-likelihood  estimate  of  the  covariance  has  been  utilized 
to  form  an  adaptive  detection  algorithm  [29],  showing  a  dramatic  increase  in  the  probability  of 
detection  compared  to  the  unconstrained  case.  This  estimator  is  useful  for  uniform  linear  arrays; 
however,  other  array  geometries  do  not  lead  to  the  Toeplitz  structure. 

The  focus  of  this  report  is  to  improve  detector  and  beamformer  performance  for  arbitrary  but 
known  arrays  and  for  an  arbitrary  number  of  interference  sources.  The  results  of  this  report  are 
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not  restricted  to  uniform  linear  arrays  with  the  corresponding  Toeplitz  covariance  matrix  structure. 
Arrays  with  arbitrary  spatial  geometries  can  also  lead  to  structured  covariance  matrices,  although 
the  stn'cture  is  not  as  obvious  as  the  Toeplitz  structure. 

This  structure  will  be  discussed  further  in  Chapter  5.  In  the  next  chapter  the  CFAR  adaptive 
matched-filter  detector  will  be  discussed. 
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4.  AN  ADAPTIVE  MATCHED-FILTER  DETECTOR 


In  this  chapter,  an  adaptive  detector  is  derived  and  the  performance  analyzed.  This  detector 
utilizes  the  Gaussian  density  function  and  the  structured  mean  but  does  not  utilize  any  information 
about  the  covariance  structure.  In  order  to  form  an  unstructured  estimate  of  the  covariance, 
it  is  assumed  that  there  are  K  independent  signal-free  data  vectors  that  are  used  to  form  an 
unconstrained  estimate  of  the  covariance  matrix.  This  estimate  is  then  substituted  into  the  known 
covariance  generalized  likelihood  ratio  test.  The  resulting  test  statistic  has  the  form  of  a  normalized 
colored-noise  matched  filter  and  is  a  CFAR  detector. 

The  purpose  in  investigating  this  detector  is  that  the  resulting  test  statistic  is  simpler  than 
that  derived  through  the  full  unconstrained  generalized  likelihood  ratio  procedure.  For  real-time 
applications,  the  difference  in  complexity  may  be  significant.  Adaptive  detectors  using  constrained 
covariance  estimates  can  be  derived  by  the  same  method  used  to  derive  the  detector  presented  in 
this  chapter.  This  detector  and  the  unconstrained  generalized  likelihood  ratio  detector  provide  a 
performance  base  to  which  the  performance  of  constrained  adaptive  detectors  can  be  compared. 

Equations  describing  the  performance  of  this  detector  are  derived  for  signals  on  boresight 
(i.e.,  in  the  d  direction)  as  well  as  for  signals  which  are  not  matched  to  boresight.  The  performance 
of  this  detector  will  be  compared  in  detail  to  other  detection  methods  in  Chapter  9. 

4.1  Derivation  of  the  Test  Statistic 

The  signal  model  is  a  slight  variation  of  that  given  in  Chapter  2.  For  notational  simplicity 
the  primary  data  vector  z  will  be  considered  separately  from  the  secondary  data,  z  is  assumed 
to  be  a  complex  7V-length  Gaussian  random  vector  with  mean  0  under  hypothesis  Hq ,  mean  6d 
under  hypothesis  H j,  and  positive-definite  covariance  matrix  R.  K  >  N  additional  independent 
data  vectors  are  arranged  as  the  columns  of  the  data  matrix  Z.  These  secondary  data  vectors 
are  assumed  to  have  mean  0  and  covariance  R.  The  restriction  on  the  number  of  secondary  data 
vectors  is  made  so  that  the  maximum-likelihood  estimate  of  the  covariance  matrix  will  exist  almost 
surely. 

The  procedure  used  to  derive  the  test  statistic  is  to  assume  that  the  covariance  is  known 
and  then  to  write  the  generalized  likelihood  ratio  test  maximizing  over  the  unknown  parameter 
b.  The  resulting  test  statistic  is  the  output  power  of  the  standard  colored-noise  matched  filter. 
The  maximum-likelihood  estimate  of  the  covariance  based  on  the  secondary  data  alone  is  then 
substituted  into  this  test. 

The  derivation  is  begun  by  writing  the  generalized  likelihood  ratio  test 
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Substituting  the  complex  multivariate  Gaussian  density  functions  and  canceling  common 
terms  yields 


Now  the  logarithm  can  be  taken,  and  the  result  simplified  to 

log(A)  =  29?(6*dtR~1z)  -  Ih^R^d 

Maximizing  this  with  respect  to  the  unknown  complex  amplitude  b  yields 

-  d'R~lz 
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Substituting  Equation  (69)  into  Equation  (68)  and  simplifying  produces  the  test 
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with  a  =  log(q).  This  test  statistic  is  proportional  to  the  squared  magnitude  of  the  output  of  the 
colored-noise  linear  matched  filter  because  the  term  in  the  denominator  is  a  constant  when  the  true 
covariance  is  known. 

If  the  noise  covariance  matrix  were  known,  then  the  detector  described  by  Equation  (70) 
would  be  used.  In  general,  the  covariance  matrix  is  unknown  and  must  be  accounted  for  by  using 
adaptive  techniques.  The  generalized  likelihood  ratio  test  (GLRT)  provides  one  such  adaptive 
approach.  We  propose  to  account  for  not  knowing  the  true  covariance  by  the  ad  hoc  procedure  of 
substituting  the  maximum- likelihood  estimate  based  on  the  secondary  data. 

R  =  yrZZ1  (71) 

A 

Reed,  Mallett,  and  Brennan  (RMB)  used  a  similar  approach  in  their  maximum  signal- to- noise 
formulation  of  the  detection  problem.  The  test  form  is  then 
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We  call  this  test  the  adaptive  matched-filter  (AMF)  test.  This  test  statistic  has  the  RMB  test 
statistic  as  the  numerator,  with  a  normalization  that  is  the  same  as  that  which  would  be  provided 
by  the  unconstrained  GLRT  for  a  large  number  of  secondary  samples.  This  normalization  will 
provide  the  desired  CFAR  behavior  and  is  a  natural  normalization  factor  to  use  for  this  purpose. 

The  AMF  test  may  also  be  derived  by  other  methods.  This  test  statistic  also  results  from  a 
type  of  cell-averaging  CFAR  where  the  cell  average  is  made  from  the  outputs  of  an  RMB  adaptive 
beamformer  [16]. 

The  GLRT  uses  all  the  data  (primary  and  secondary)  in  the  likelihood  maximization  under 
each  hypothesis.  The  AMF  test  makes  no  use  of  the  primary  vector  to  estimate  the  covariance; 
therefore,  poorer  detection  performance  might  be  expected.  In  the  following  sections,  the  per¬ 
formance  loss  is  shown  to  be  small  and  that,  in  certain  situations,  the  AMF  test  will  actually 
outperform  the  GLRT. 

4.2  CFAR  Behavior 

We  now  show  that  the  density  of  the  AMF  test  statistic  does  not  depend  on  the  true  covariance 
matrix  under  Hq.  and  thus,  it  gives  a  constant  false  alarm  rate  test. 

Let  u  =  R~id.  and  y  =  R~^z.  Then  the  test  can  be  written 


lu+RiR^Riy!2  I^R-’-y!2 

- — - - — —  =  - r -  £  Q  ,  (7o) 

u'RjR_1R5 u  u*R  ;u  ^ 

where  R  =  R~^RR~i.  R  is  subject  to  the  complex  Wishart  distribution  with  parameters  K ,  N. 
and  I,  which  will  be  denoted  CW(A',.Y;I)  [5], 

Now  a  unitary  transform  is  defined  that  rotates  the  whitened  signal  vector  into  the  first 
elementary  vector. 


de  =  L^u.  e  =  [1.0 . 0[t  (74) 

The  first  column  of  U  is  the  whitened  signal  vector  u,  and  the  other  N  -  1  columns  form  an 
arbitrary  orthonormal  basis  for  the  orthogonal  complement  of  the  subspace  spanned  by  u.  The 
test  then  becomes 
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where  x  =  Ufy  and  S  =  TJtRU.  Then  x  is  distributed  N(0,I)  under  Hq,  and  S  is  distributed 
CW(K,  N;l).  The  actual  covariance  does  not  appear  in  this  equation  or  in  the  underlying  density 
functions,  and  thus,  this  is  a  CFAR  test.  This  test  is  independent  of  both  the  scale  and  the 
structure  of  the  true  covariance,  in  contrast  to  the  simple  unknown  level  CFAR  characteristic  of 
many  common  CFAR  detectors. 

4.3  Generalization  of  the  Signal  Model 

The  equations  to  determine  the  performance  of  this  detector  will  be  derived  for  a  general 
signal  case  where  the  signal  may  or  may  not  be  in  alignment  with  the  look  direction.  In  our  model 
the  signal  is  assumed  to  lie  along  some  general  direction  vector  p;  hence,  the  signal  is  now  normally 
distributed  N(0,  R)  on  Hq  and  N(  bp,  R)  on  H\.  The  steering  vector  of  the  array  is  assumed  to 
be  q. 

The  direction  vectors  may  be  normalized  so  that 

pfp  =  q+q  =  i  ,  (76) 

and  the  following  definitions  are  made: 

A2q  =  (qtR^q)  ,  (77) 

and 

A2p  =  (ptR_1p)  (78) 

Summarizing  Kelly  [30],  these  terms  may  be  used  to  describe  the  signal-to-noise  ratio.  That 
is,  the  maximum  signal-to-noise  ratio  is 

SNRqq  =  \b\2Aq  ,  (79) 


which  is  attained  when  the  signal  lies  along  the  axis  for  which  the  detector  is  steered.  There  is  a 
signal-to-noise  ratio  loss  when  the  signal  does  not  lie  in  the  steering  direction.  The  signal-to-noise 
ratio  which  results  when  the  array  is  steered  in  the  direction  corresponding  to  q  is 
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q^R  *p  can  be  interpreted  as  an  inner  product  of  p  and  q,  and  the  definition  that 


cos  deJ<?  = 


ApAq 


(81) 


can  be  made  [21]. 


Then  cos  6  may  be  used  to  relate  the  signal-to-noise  ratio  to  the  maximum  signal-to-noise  ratio 


as 


a  =  SNRqp  =  SNRpp  cos2  6  (82) 

STVflqp  can  be  thought  of  as  the  signal-tonoise  ratio  in  the  subspace  spanned  by  the  adapted 
steering  direction,  and  likewise, 

c  =  SNRpp  sin2  8  (83) 

can  be  viewed  as  the  signal-to-noise  ratio  in  the  orthogonal  subspace. 

4.4  Performance  Evaluation 

4.4.1  Derivation  of  Test  Performance 

The  analysis  of  this  detector  is  similar  to  Kelly's  analysis  [14]  and  uses  the  same  notation.  Ap¬ 
propriate  whitening  and  unitary  transforms  as  defined  in  Section  4.2  are  performed  to  reformulate 
the  AMF  test  in  the  statistically  equivalent  form 
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In  Equation  (75)  z  has  been  redefined  in  this  equation  to  be  the  whitened  rotated  primary  data 
vector  x. 


3This  definition  is  equivalent  to  the  definition  of  available  signal-to-noise  ratio  a  in  Equation  (44) 
when  p  =  q. 
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The  steps  made  to  form  this  representation  are  identical  to  those  used  to  show  that  the 
test  statistic  is  independent  of  the  underlying  covariance  matrix.  Here,  because  of  the  generaliza¬ 
tion  of  the  signal  model,  the  actual  mean  direction  vector  may  not  have  been  transformed  to  the 
first  elementary  vector.  With  this  in  mind,  z  is  now  normally  distributed  N(0, 1)  under  Ho  and 
N(Mpf,  I)  under  Hi.  The  transformed  covariance  estimate  S  has  the  complex  Wishart  distribution 
CW(K,  AT;  I),  and  the  signal  direction  vector  is  given  by 


f  =  -J-UfR"ip 


(85) 


U  is  the  unitary  matrix  required  to  rotate  the  whitened  direction  vector  to  the  first  elementary 
vector  e  =  [1,0,...,  0]*.  Following  Kelly’s  method  [14],  the  vector  f  is  decomposed  into  components 
parallel  and  orthogonal  to  e,  decomposing  S  as  well: 
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The  test  statistic  can  now  be  simplified  to  the  form 
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and  thus,  the  test  may  be  expressed 
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where  v  is  now  normally  distributed  N(0, 1)  or  N(v/pi4p6cos0, 1). 
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The  GLRT  is  expressed  in  a  similar  form  without  the  p  factor  in  the  threshold.  T  is  an 
independent  random  variable  which  is  distributed  chi-squared  (x2)  with  L  complex  degrees  of 
freedom.  Equation  (91)  has  the  form  of  a  scalar  CFAR  test,  in  which  the  threshold  is  multiplied 
by  the  random  loss  factor  p. 

The  density  of  the  loss  factor  p  has  been  derived  [30],  and  it  is  given  by  the  formula 


f(p)  =  e~cp 


(N  +  L-  1)! 
(iV  +  L  —  1  +  m)! 


cmfj(p;L  +  l.N  +  m-1) 


(92) 


where  c  is  given  in  Equation  (83),  and  we  have  defined  L  =  K  4-  1  —  N .  The  central  beta  density 
function  is  defined  by 


/j(i;  n,  m)  =  (”  ^  xn~\l  -  x)^1 

( n  -  l)'(m  —  1): 


(93) 


An  alternative  form  for  f(p)  may  be  derived  by  expressing  this  function  in  terms  of  the 
confluent  hypergeometric  function  and  using  Rummer’s  first  transformation  [31]  to  yield  [30] 


x  cm 


f(p)  =  e  c  £  L-v  +  m  -  1) 

‘  771! 


(94) 


m=0 


4.4.2  Evaluation  of  the  Probability  of  False  Alarm 

The  probability  of  false  alarm  (PFA)  for  the  AMF  test  is  calculated  when  the  signal  mean  is 
equal  to  0:  consequently,  the  orthogonal  SNR  term  c  is  0,  and  the  density  functions  for  p  reduce  to 
the  central  beta  density  function.  The  probability  of  false  alarm  will  have  the  same  form  as  that 
of  the  GLRT  [15,14]  except  for  the  presence  of  the  factor  p  in  the  threshold.  As  shown  [21],  the 
probability  of  false  alarm  for  the  GLRT  is  given  by 


PFAglrt  = 


1 


(1  -i-  a)' 


(95) 


where  a  =  7/(1  -  7),  and  7  is  the  threshold  term  of  Equation  (63).  To  determine  the  false  alarm 
probability  for  the  AMF  test,  the  term  a  can  be  replaced  with  pa,  and  the  expectation  with  respect 
to  the  loss  factor  p  can  be  taken  to  yield 


PFA 


AMF 


=  [  PF Aamf\p  f{p)dp  -  [ 
Jo  Jo 


'Mp-.L  +  LN-  1) 
(1  +  ap)L 


dp 


(96) 
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This  has  been  evaluated  through  numerical  integration,  and  also  by  means  of  an  expansion  into 
ait  infinite  series,  integration  term  by  term,  and  derivation  of  truncation  bounds,  with  the  same 
results.  Iterative  procedures  based  on  bisection  and  Newton's  method  have  been  used  to  find  a 
when  a  particular  PFA  is  specified. 

4.4.3  Evaluation  of  the  Probability  of  Detection 

The  conditional  probability  of  detection  for  the  AMF  test,  given  p,  may  be  expressed  in  a 
finite  sum  expression  as  [32] 


PDamf\p  =  1  - 


(1  +  ap)L 


to 

m=  I 


MmG„ 


ap 


1  +  ap 


(97) 


where  a  is  the  SNR  component  defined  in  Equation  (82)  and  where  Gm  is  the  incomplete  Gamma 
function 


m—  1  k 

Gm(y)  = 

k= 0  K' 


(98) 


Unlike  the  GLRT,  the  AMF  test  includes  the  loss  factor  in  the  threshold  as  well  as  in  the  mean. 
The  expectation  of  the  conditional  probability  of  detection  with  respect  to  the  random  variable  p 
must  be  taken  to  evaluate  the  unconditional  probability  of  detection.  The  probability  of  detection 
can  be  written  as 


pD^=i-j!(T^ts™y°prGm(^,{p)dp  ■  <99) 

This  equation  has  been  computed  through  numerical  integration  using  the  finite  sum  form 
of  the  density  function  of  Equation  (92).  Additionally,  this  equation  has  been  evaluated  through 
the  use  of  the  infinite  series  form  of  the  density  function  by  integrating  term  by  term  to  express 
the  probability  of  detection  as  a  series  expression  containing  two  finite  and  three  infinite  series. 
Bounds  for  the  three  infinite  series  were  obtained  using  methods  similar  to  that  of  Shnidman  [33]. 
The  results  using  this  method  were  then  used  to  verify  the  results  of  the  numerical  integration. 

In  order  to  evaluate  the  probability  of  detection  numerically,  a  single  routine  has  been  written 
to  evaluate  the  probability  of  detection  of  a  scalar  CFAR  detector  [32]: 


PDcfar{a,  a,  L)  =  1  - 


(1 


'  m  =  i 


(100) 
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The  numerical  integration  technique  of  finding  the  unconditional  probability  of  detection  consists 
of  repeatedly  calling  the  routine  with  a  and  a  replaced  by  ap  and  ap.  weighting  the  result  by  /j(p), 
and  summing.  When  this  procedure  is  performed  for  the  AMF  test,  the  equation  implemented  is 

7-1  1 

PD  AMF  =  ApJ2  PDcfar(iApa ,  iApa.  L)fp{iAp),  Ap  =  -  ,  (101) 

i—0 

with  a  defined  in  Equation  (82)  as  the  SNR  component  parallel  to  the  direction  vector.  I  is  chosen 
to  yield  a  suitably  small  error  by  successively  doubling  the  number  of  terms  until  the  probability 
of  detection  varies  less  than  some  e.  It  typically  requires  between  1024  and  8192  integration  terms 
to  achieve  an  error  bound  of  0.00001. 

The  corresponding  equation  for  the  GLRT  is  numerically  integrated  in  the  same  manner,  with 
the  equation  implemented  being 

i-i 

PDgirt  =  Ap  ]T  PDcfar{a ,  iApa.  L)fp(iAp)  .  (102) 

t=0 

If  a  in  Equations  ( 101 )  or  ( 102)  is  now  replaced  by  7a.  where  7  is  a  random  loss,  the  probability 
of  detection  for  the  Swerling  target  fluctuation  models  [3]  may  be  found.  For  these  cases,  7  is  subject 
to  a  \2  distribution  with  the  number  of  degrees  of  freedom  dependent  upon  the  Swerling  model 
chosen.  For  the  Swerling  I  model  [3j.  there  will  be  only  one  complex  degree  of  freedom,  and  the 
density  function  for  7  is 

f,  =  e"  •  (103) 
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PDclrt\P  = 


1  +  ap 
1  +  a  +  ap 


L 


(106) 


In  these  cases,  no  signal  mismatch  is  assumed.  For  comparison,  the  probability  of  detection  for  the 
Swerling  I  known  covariance  matched  filter  is  [32] 


PDmf  =  e_Q/(1+a)  .  (107) 

In  Figure  7,  it  can  be  seen  that  the  AMF  detector  requires  a  slightly  higher  signal-to-noise  ratio 
to  achieve  the  same  probability  of  detection  as  the  unconstrained  GLRT  detector  for  lower  signal-to- 
noise  ratios,  with  a  crossover  for  the  probability  of  detection  of  0.9.  The  known  covariance  detector 
performs  better  than  either  of  the  two  adaptive  approaches.  These  detectors  may  be  utilized  if 
the  increase  in  signal-to-noise  ratio  required  to  achieve  the  performance  of  the  known  covariance 
detector  is  acceptable.  When  the  quantity  of  data  is  limited,  the  unconstrained  detectors  will  not 
exist  since  the  unconstrained  maximum-likelihood  estimate  of  the  covariance  matrix  requires  that 
K>N .  For  K  =  N ,  tens  of  decibels  more  signal-to-noise  ratio  may  be  required  to  provide  the  same 
detection  performance  as  the  known  covariance  detector,  as  is  shown  in  the  first  of  the  plots.  For 
a  probability  of  detection  of  0.8,  the  GLRT  and  the  AMF  detectors  require  approximately  11  and 
12  dB  greater  signal-to-noise  ratio  for  N  =  4,  K  =  4.  This  loss  in  performance  is  unacceptable  in 
most  applications  and  motivates  us  to  use  our  knowledge  of  the  structure  of  the  covariance  matrix 
to  increase  the  adaptive  detector  performance. 

The  basic  idea  found  in  the  derivation  of  the  AMF  detector  will  be  utilized  again  when  a 
constrained  estimate  of  the  covariance  matrix  is  substituted  into  the  known  covariance  matched- 
filter  detector.  In  Chapter  9,  the  performance  of  detectors  based  on  constrained  and  unconstrained 
covariance  estimates  will  be  compared. 
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Figure  7.  Probability  of  Detection  versus  SSR.  Known  covariance  detector  compared  to 
unconstrained  GLRT  and  AMF  detectors. 
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5.  COVARIANCE  STRUCTURE 


5.1  Introduction 

Th°  complete  signal  model  derived  in  Chapter  2  is  not  commonly  used  in  the  derivation  of 
adaptive  signal-processing  algorithms.  Estimates  of  the  covariance  matrix,  used  to  form  the  test 
statistic  or  weight  vector,  are  typically  not  restricted  to  covariance  matrices  that  are  feasible  with 
a  particular  array  or  noise  environment.  In  subsequent  chapters,  the  complete  signal  model  will 
be  used  in  the  derivation  of  the  adaptive  beamforming  and  detection  algorithms.  In  this  chapter, 
covariance  matrices  that  are  consistent  with  the  signal  model  are  characterized,  and  methods  of 
representing  the  matrices  are  presented.  The  parameters  provided  by  the  representations  enable 
the  use  of  the  maximum-likelihood  procedure  to  estimate  covariance  matrices  that  are  constrained 
as  a  result  of  the  data  model. 

The  noise  at  the  output  of  the  sensor  is  assumed  to  be  the  sum  of  two  independent  terms. 
The  first  term  is  receiver  noise,  which  is  assumed  to  be  independent  from  sensor  to  sensor.  The 
contribution  of  this  noise  to  the  output  covariance  matrix  is  then  a  diagonal  matrix  and  may  be 
known.  The  effect  of  this  matrix  on  the  structure  of  the  sensor  output  noise  covariance  matrix  will 
be  investigated  later  in  this  chapter. 

The  second  noise  term  is  due  to  the  radiation  received  by  the  sensors.  The  responses  of  the 
sensors  will  not.  in  general,  be  mutually  independent  since  the  cross  covariance  between  the  received 
signals  is  assumed  to  satisfy  the  signal  model  found  in  Equation  (32)  and  developed  in  Chapter  2. 
The  contribution  of  this  noise  term  to  the  output  covariance  matrix  will  lie  within  a  class  of  matrices 
ns.  This  class  is  characterized  by  an  integral  expression  involving  a  non-negative  spatial-temporal 
power-spectral-density  function  of  Equation  (32).  Non-negativity  of  the  spatial-temporal  power- 
spectral-density  function  and  linearity  of  the  integral  expression  imply  that  the  covariance  matrix 
will  exhibit  a  structure  that  is  dependent  upon  the  spatial  locations  of  sensors,  the  spatial  and 
temporal  response  of  the  sensors,  and  the  transmission  medium.  The  actual  covariance  matrix  will 
be  dependent  upon  the  intensity  and  location  of  the  spatially  distributed  noise  sources. 

The  class  of  possible  covariance  matrices  will  be  examined  from  geometric  and  algebraic 
viewpoints.  This  provides  methods  of  representing  the  covariance  matrices  by  the  use  of  a  finite 
number  of  terms  and  provides  the  base  upon  which  adaptive  signal-processing  algorithms  will  be 
built. 

5.2  Characterization  of  the  Constraint  Space 

5.2.1  Constraint  Space  Definition 

The  covariance  matrix  has  been  shown  to  result  from  an  integral  expression  involving  a 
non-negative  spatial-temporal  power-spectral  density  along  with  functions  that  are  defined  by  the 
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sensor  responses  and  the  geometry  of  Equation  (32).  The  set  of  all  matrices  that  satisfy  the  spatial 
constraint  will  be  denoted  by 

fts  =  {R  :  R  =  fe,#  Svv(u0, 9 , 0)d(u>o,  9 ,  0)df{uo, 9, 0)  cos  9d9d0, 

Svv{u>0,9,<t>)>  0}  .  (108) 

We  w’ould  now  like  to  characterize  the  set  of  covariance  matrices  that  are  members  of  Qs. 

Covariance  matrices  that  are  members  of  fi,  belong  to  a  larger  vector  space  whose  elements 
are  general  N  by  N  Hermitian  matrices.  For  Hermitian  matrices  Ri  and  R2  and  real  scalars  a  and 
3.  qRi  +  ^3R2  is  a  Hermitian  matrix,  and  this  can  be  used  to  define  a  vector  space.4 

An  inner  product, 

<  Rj,R2  >=  tr(RiR2)  ,  (109) 

can  be  defined  on  this  space,  and  for  completeness,  the  obvious  norm  is 

j|R|i2  =  tr(RR)  (110) 

It  will  be  convenient  to  write  the  elements  of  the  covariance  matrix  R  as  an  N2  component 
real  vector  [35], 

r  =  [rj  . . .  rV2jr  .  (Ill) 

Let 

r,  =  Ra  1  <i<  N  ,  (112) 

and  let 

CjV+i  +  jry+2  =  V  2R\2  ■  (113) 


4  Although  I  have  used  R  to  refer  to  the  elements  of  this  space,  the  members  of  this  space  include 
all  Hermitian  matrices,  not  just  non-negative  definite  matrices. 
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By  continuing  this  process  along  the  upper  minor  diagonals  of  R  in  a  generalized  store-by-diagonal  [36] 
ordering,  then 


r  =  \Rn  •  -  •  Rnn  V2R(R12)  y/2 $(RU)  ■  ■  ■  y/2 R(fl1JV)  y/2Qs{R1N)]T  .  (114) 

The  precise  ordering  of  the  elements  of  r  is  unimportant  to  the  results  that  utilize  this  notation. 
Because  of  Hermitian  symmetry,  R  is  fully  specified  by  writing  its  upper  triangle  in  this  form. 

With  this  definition  of  r,  there  is  an  equivalence  between  the  familiar  vector  inner  product 
and  the  inner  product  given  by  Equation  (109). 


<  RltR2  >=  tr(RiR2)  =  rj'r2 


(115) 


5.2.2  Convex  Cone 

The  space  Q,  of  covariance  matrices  is  a  convex  cone.  This  follows  directly  from  the  linearity 
of  the  integral  and  the  positivity  of  the  spatial-temporal  power-spectral  density.  To  see  this,  let  5j 
and  S2  be  valid  distributions  of  spatial  energy  with  covariance  matrices  Rj  and  R2,  respectively, 
Rj ,  R2  €  Then.  53  =  aSi  +  3S2  is  a  valid  spatial  energy  distribution  for  real  o.  3  >  0.  The 
linearity  of  the  covariance  integral  then  implies  that  R3  =  aR]  -I-  3R-2  is  a  valid  covariance  matrix, 
so  R3  £  f is. 

The  properties  of  convex  sets  can  be  used  to  parameterize  the  covariance  matrices  that  are 
in  our  constraint  space. 

5.2.3  Dimension  of  Space 

The  space  of  possible  covariance  matrices,  ft,,  forms  a  subset  within  the  set  of  Hermitian 
matrices.  A  characteristic  of  the  subset,  dependent  upon  the  sensor  responses  and  array  geometry, 
is  the  dimension  (L)  of  the  space.  The  dimension  L  of  Q,  is  the  maximum  number  of  elements  of 
Q,  that  is  linearly  independent  for  real  coefficients.  That  is,  for  I  >  L  there  will  exist  a  solution  to 

/ 

]Ta,R,=0  (116) 

t=i 

for  real  a,  not  all  zero. 

We  can  readily  see  that  the  maximum  dimension  of  ft,  will  be  N2.  We  will  make  the  definition 
that  the  covariance  matrices  for  an  array  exhibit  structure  if  the  dimension  of  Qa  is  less  than  N2; 
n,  is  a  convex  cone  within  a  proper  subspace  of  the  N2  real  vector  space.  It  will  be  assumed  that 
the  array  is  such  that  the  covariance  matrices  will  be  structured. 
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5.2.4  Examples 


The  minimum  dimension  of  fi3  for  unknown  covariance  matrices  occurs  when  the  noise  envi¬ 
ronment  is  known  except  for  a  scale  factor.  The  dimension  of  this  space  is  1.  A  CFAR  detector 
based  on  this  minimum  constraint  space  is  derived  and  analyzed  in  Appendix  A.  Closed-form 
expressions  result  for  the  test  statistic  and  the  equations  that  describe  the  performance  of  this 
detector.  Because  the  model  for  the  unknown  covariance  matrix  in  this  case  is  the  most  restrictive 
constraint-cone  possible,  the  performance  of  this  test  provides  a  convenient  reference  to  which  the 
performance  of  other  constrained  detectors  can  be  compared. 

It  will  be  instructive  to  consider  examples  of  what  the  dimension  would  be  for  other  array 
geometries.  The  maximum  dimension  for  sensors  with  identical  spatial  response  is  N2  -  N  + 
1.  Covariance  matrices  for  this  array  will  have  equal  diagonal  elements  but  otherwise  may  have 
considerable  freedom  for  the  values  of  the  off-diagonal  terms.  The  maximum  dimension  occurs 
for  sensor  spacings  such  that  there  are  no  repeated  off-diagonal  covariance  “lags.”  This  would  be 
quite  easy  to  achieve  since  sensor  spacings  are  continuous  quantities.  Sensors  can  be  spaced  so  that 
the  three-dimensional  differences  in  the  array  positions  do  not  repeat.  This  set  of  array  position 
differences  is  known  as  the  co-array  [37]. 

A  linear  array  whose  array  spacings  are  a  multiple  of  a  common  factor  and  whose  covariance 
matrix  will  contain  no  repeated  or  missing  off-diagonal  lags  is  known  as  a  zero  redundancy  array.  A 
zero  redundancy  array  cannot  always  be  formed  from  N  array  elements.  An  array  with  the  minimum 
number  of  repeated  lags  is  known  as  a  minimum  redundancy  array.  Pillai  [38]  has  tabulated  the 
array  spacings  for  zero  and  minimum  redundancy  arrays.  Some  of  the  simulation  results  shown  in 
Chapter  9  will  be  for  a  minimum  redundancy  array. 

For  uniform  linear  arrays,  we  have  seen  in  Section  2.7.2  that  the  underlying  covariance  matrix 
will  exhibit  the  Toeplitz  structure.  Here,  the  dimension  is  2 N  -  1.  A  detection  method  based  on 
restricting  covariance  matrices  to  this  structure  has  been  investigated  by  Fuhrmann  [29]. 

If  the  eigenvectors  of  the  covariance  matrix  are  known,  then  there  is  a  unitary  transformation 
that  can  be  applied  to  the  data  so  that  the  covariance  matrix  is  diagonal.  This  situation  occurs 
when  the  spatial  noise  field  is  assumed  to  be  periodic,  with  the  period  given  by  the  array  aperture 
for  a  uniform  linear  array.  In  this  instance,  the  dimension  of  the  space  will  be  N.  Fuhrmann  [39] 
has  investigated  detection  algorithms  based  on  this  constraint  space  as  well. 

5.3  Representations 

In  order  to  develop  adaptive  signal-processing  algorithms  that  make  use  of  the  constraints 
induced  by  the  signal  model,  methods  of  enforcing  the  constraints  while  performing  the  parameter 
estimation  must  be  developed.  The  methods  that  are  developed  here  are  based  on  two  approaches. 
The  first  approach  is  to  form  a  representation  for  the  entire  space  of  possible  covariance  matrices. 
Approximating  the  integral  equation  that  defines  the  constraint  space  by  a  finite  sum  is  the  method 
utilized  to  generate  this  representation  for  the  constraint  space. 
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The  second  approach  is  through  the  representation  of  a  covariance  matrix  as  a  member  of  a 
convex  set.  The  characterization  of  the  constraint  space  as  a  convex  set  allows  certain  representa¬ 
tions  for  covariance  matrices  that  are  members  of  the  set. 

Both  of  the  approaches  develop  representations  parameterizing  the  covariance  matrices  with  a 
finite  number  of  terms.  The  finite  representations  allow  the  practical  implementation  of  covariance 
estimators  with  estimates  constrained  to  be  within  ft,  of  Equation  (108). 

5.3.1  Hyperplane  Characterizations 

Before  further  discussion  of  the  representations,  the  algebraic  and  geometric  idea  of  a  hyper- 
plane  will  be  introduced.  A  hyperplane  H  in  a  linear  vector  space  ft  can  be  defined  as  [40,41  j5 

H  =  {R:  tr(RG)  =  0,  G  #  0,  G  =  G+;  G,R  6  ft}  (117) 

The  hvperplane  that  results  from  this  definition  is  a  maximal  proper  subspace  of  ft.  This  equa¬ 
tion  defines  an  algebraic  subspace,  which  is  orthogonal  to  G.  The  space  of  possible  constrained 
covariance  matrices  is  a  reduced  dimensional  subset  of  all  possible  covariance  matrices  and  is  then 
a  subset  of  a  hyperplane  of  the  larger  dimensional  space. 

Hvperplanes  are  of  interest  because  the  boundaries  of  cones  formed  from  a  finite  number  of 
terms  are  hyperplanes  or  intersections  of  hvperplanes.  The  inner  product  of  any  element  on  the 
hyperplane  boundary  of  the  constraint  set  with  at  least  one  of  the  elements  of  the  vector  space  will 
be  zero. 

Additionally,  there  are  some  representations  for  covariance  matrices  that  can  be  viewed  as 
hyperplane  representations. 

5.3.2  Finite  Sum  Representation 

The  first  representation  that  will  be  used  to  parameterize  the  constraint  space  is  to  approxi¬ 
mate  the  covariance  integral  found  in  Equation  (28)  by  an  A/-term  finite  sum6 

At 

ft,  =  {R  :  R  =  "^2  Om,Qm)d{u;o,6rn.0m)df(^o,^m,Om)}  i  (118) 

m  =  1 


5More  generally.  tr(RG)  =  c,  an  arbitrary  constant.  Our  attention  is  placed  on  hvperplanes  that 
pass  through  the  origin. 

6The  differential  area  or  volume  is  a  positive  quantity,  which  can  be  incorporated  into  the  scaling 
of  the  direction  vectors  or  the  spatial  spectrum. 


49 


where  the  number  of  terms  M  and  the  locations  of  the  spatial-temporal  spectral  points  (0m><Pm) 
have  been  chosen  to  produce  a  suitable  approximation  to  the  integral  of  Equation  (28).  One  such 
distribution  of  points  is  spatially  uniform  in  (#,  q).  The  effect  of  this  is  to  approximate  the  convex 
cone  of  the  constraint  set  by  an  inscribed  polyhedral  cone  [42,37].  This  representation  will  be  used 
in  Chapter  6. 

The  outer  products  of  the  direction  vectors  are  basis  functions  in  the  sense  that  any  covariance 
matrix  interior  to  the  polyhedral  cone  can  be  formed  by  the  weighted  sum  of  these  basis  functions. 
These  functions  are  not,  in  general,  independent  and  may  not  form  an  algebraic  basis  for  the 
covariance  matrices  that  are  in  the  constraint  set.  The  definition  of  the  cone  requires  that  the 
weights  for  the  basis  functions  be  non-negative;  the  algebraic  basis  does  not  make  this  restriction. 

If  the  bases  are  fixed,  then  not  all  matrices  that  are  members  of  the  constraint  class  can  be 
represented  in  this  manner  for  finite  M.  If  the  covariance  matrix  contains  the  contribution  from 
a  single  large  interference  source,  then  it  may  be  that  there  is  not  a  representation  of  it  in  terms 
of  the  discrete  spectrum.  For  a  finite  number  of  spectral  points,  the  set  of  covariance  matrices 
that  cannot  be  represented  is  non-empty.  If  the  number  of  points  in  our  approximation  M  is 
increased  and  an  upper  bound  is  placed  on  /  5(e)  de,  then  the  continuity  of  d(u,'0. 6,  <t>)  allows  the 
representation  to  get  arbitrarily  close7  to  any  matrix  within  f ls.  An  analogy  to  this  would  be  to 
approximate  a  disk  by  the  interior  of  an  inscribed  polygon.  As  the  number  of  sides  on  the  polygon 
is  increased,  the  area  of  the  disk  that  is  not  contained  within  the  polygon  can  be  made  arbitrarily 
small.  Figure  8  illustrates  this  effect  for  the  covariance  matrix  with  two  identical  sensors  spaced 
one-half  wavelength  apart.  If  the  covariance  matrices  possible  when  R\i  =  1  are  examined,  then 
i?i2  will  be  constrained  bv  the  requirement  that  the  covariance  matrix  be  non-negative  definite. 
The  possible  locations  of  Ru  will  be  constrained  to  lie  within  a  disk  such  that  |i?i2|  <  1-  The 
approximation  to  this  constraint  space  for  four  terms  in  the  finite  sum  is  shown  by  the  inscribed 
square.  This  example  will  be  investigated  in  more  detail  in  Section  5.5. 

It  will  be  assumed  that  a  large  number  of  terms  will  be  required  to  parameterize  covariance 
matrices  that  are  members  of  Spatial  receiving  arrays  often  have  high-powered  discrete  inter¬ 
ference  sources;  the  covariance  matrices  that  result  are  poorly  conditioned.  If  the  spatial  spectrum 
is  known  to  be  “smooth,”  then  this  information  could  be  used  to  reduce  the  number  of  terms  in 
the  summation. 

5.3.3  Caratheodory  and  Maximal  Representations 

The  large  number  of  terms  used  in  the  previous  section  to  parameterize  the  constraint  space  is 
not  necessary  to  represent  a  single  covariance  matrix  that  is  a  member  of  this  constraint  space.  The 
spectral  weights  of  the  previous  section  are  not  generally  unique  due  to  the  linear  dependence  of  the 


7The  norm  of  the  difference  between  any  covariance  matrix  within  fl,  and  the  nearest  covariance 
matrix  within  the  polyhedral  cone  can  be  made  arbitrarily  small. 
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Figure  8.  Constraint  regions  and  approximations. 


basis  functions.  This  algebraic  property  of  the  basis  functions  can  be  used  to  find  a  representation 
of  the  covariance  matrix  requiring  fewer  terms.  The  number  of  terms  that  will  be  used  in  this 
representation  is  dependent  upon  the  dimensionality  of  the  constraint  space.  In  the  previous  section 
any  member  of  the  constraint  set  could  be  represented  arbitrarily  closely  by  adjusting  the  spectral 
weights.8  This  will  not  be  true  with  the  representations  introduced  here.  Rather  than  adjusting 
the  weights  of  a  large  number  of  basis  functions,  the  weights  and  the  directions  of  arrival  of  a  much 
smaller  number  of  bases  can  be  adjusted.  This  parameterization  of  the  members  of  the  constraint  set 
can  be  used  to  estimate  structured  covariance  matrices.  The  intent  here  is  to  estimate  a  structured 
covariance  matrix  with  a  reduced  number  of  terms.  To  accomplish  this  the  polyhedral  cone  formed 
by  the  finite  basis  must  be  allowed  to  ‘'move’’  to  enclose  the  estimate. 

Theorem  5.1  Caratheodory  Representation  [40.41.43,44]  Every  element  R  in  an  L-dimensional 
space  f l  that  is  a  finite  convex  combination  of  elements  of  Q  can  be  represented  as  a  convex  combi¬ 
nation  of  L  +  1  or  fewer  elements  ofQ. 


Proof  [441. 

The  proof  of  this  theorem  is  instructive  because  it  is  one  method  of  reducing  the  large  number  of 
terms,  used  in  the  finite  sum.  to  a  smaller  number  of  terms.  Using  the  vector  representation  for 
the  elements  of  the  vector  space,  form  augmented  vectors 


r 


1 

l 


(119) 


8  As  the  number  of  bases  approaches  infinity. 
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Let  r  €  fl  be  represented  as 


M  M 

r  =  ]T<Ttrt,  ^(j,  =  1,  crt  >  0  (120) 

«=i  t=i 

If  the  number  of  terms  M  is  less  than  L  +  2,  the  proof  is  complete,  and  there  is  no  need  to  proceed 
further.  Now  let 

r '  =  XI  (121) 

i  =  l 

Since  M  >  L  +  1  then  the  r'’s  are  not  independent,  and  there  exist  ^,’s  that  are  not  all  zero  such 
that 


M 

5>.r'=0  (122) 

i  =  l 

The  first  element  of  each  vector  is  1,  so  m  =  0.  and  at  least  one  m  is  positive.  Let  a  be 
the  largest  number  such  that  a /j.t  <  at,  Vi.  a  is  finite  since  at  least  one  fi,  is  positive.  Now  let 
cr,  =  a,  -  an,.  Then 

A/  A  f  M 

=  Hcr«r<  “  QX!^‘ri  =  r  -  (123) 

1=1  i=i  t=i 

and  at  least  one  dl  =  0.  r  has  then  been  expressed  in  fewer  elements.  This  argument  can  be 
repeated  until  R  has  been  expressed  as  a  positive  linear  combination  of  at  most  L  ■+■  1  elements  of 
LI.  The  method  used  h- 1<-  to  «iuec  fhr’  number  of  terT  ~  is  known  as  the  reduction  theorem  [44]. 

This  theorem  is  easily  generalized  to  a  finite  positive  combination  of  the  elements  with 
Si=iai  =  ci  an  arbitrary  positive  constant.  The  convexity  of  the  sum  was  not  utilized  in  the 
proof.  Additionally,  when  r  is  formed  from  the  elements  of  a  convex  cone,  r,  may  be  scaled  so  that 
the  condition  5Zl=i  <?,  =  1  holds. 

If  one  of  the  terms  for  all  of  the  vector  elements  is  positive  and  non-zero,  then  there  is  no  need 
to  augment  the  vectors.  The  method  used  to  prove  the  Caratheodorv  representation  can  be  applied 
directly  to  the  vectors  rather  than  to  r',  with  the  result  that  the  vector  r  can  be  represented 
using  L  elements  rather  than  L  - 1-1. 

For  sensor  arrays  where  at  least  one  sensor  has  a  non-zero  response  throughout  space,  one 
of  the  diagonal  entries  of  the  dyads  in  Equation  (118)  will  be  positive  for  all  terms  in  the  sum. 
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A  covariance  matrix  that  is  formed  by  the  positive  weighted  sum  of  M  elements  of  Q,  can  be 
represented  as  a  combination  of  at  most  L  elements  of  Cl3. 

Representing  a  member  of  a  convex  set  in  terms  of  L  non-zero  weighted  bases  is  a  maximal 
representation  [44]  of  that  member. 

5.3.4  Hyperplane  Representation 

There  is  an  additional  Caratheodory  representation  theorem  presented  in  Grenander  and 
Szego  [45]  that  has  been  used  to  show  that  for  Toeplitz  matrices  there  is  a  representation  in  terms 
of  N— 1  ‘‘spectral’’  points  plus  a  diagonal  term  [38,46].  Thus,  any  Toeplitz  matrix  can  be  represented 
as 


■Y-l 

R  =  £  atd,d]  +  o0l  ,  (124) 

1=1 

where  the  d,  terms  are  Vandermonde  vectors,9  and  cr,  >  0.  The  spectral  points  used  in  this  repre¬ 
sentation  are  not  necessarily  the  same  as  the  spectral  Doints  introduced  in  the  finite  representation 
of  the  constraint  space  in  Section  5.3.2. 

This  is  not  the  same  Caratheodory  representation  theorem  as  that  presented  in  Section  5.3.3, 
but  it  can  be  related  if  the  representation  of  the  <7oI  term  is  investigated.  The  term  a0 1  can  be 
formed  by  the  sum  of  the  outer  products  of  N  orthogonal  Vandermonde  vectors,  the  discrete  Fourier 
transform  vectors.  Then  R  is  the  positive  weighted  sum  of  2A'  -  1  —  L  direction  vectors. 

This  representation  can  be  interpreted  for  our  problem  as  a  hyperplane  representation.  There 
will  be  at  least  one  vector  in  the  L-dimensional  vector  space  that  is  orthogonal  to  the  convex  cone 
generated  by  the  matrices  I  and  the  d,d^.  The  intersection  of  this  hvperplane  and  the  convex  cone 
contains  the  covariance  matrix  R. 

This  representation  would  appear  to  be  a  suitable  representation  for  covariance  matrices  that 
are  members  of  the  constraint  space  generated  by  a  uniform  linear  array.  Any  matrix  within  the 
constraint  space  will  have  a  representation  in  this  form;  however,  the  converse  is  not  true.  Matrices 
that  can  be  represented  in  this  form  may  not  be  members  of  the  constraint  space.  The  Vandermonde 
vectors  that  are  used  in  this  representation  may  not  be  realizable  as  direction  vectors  for  a  particular 
array  geometry  and  operating  frequency,  and  thus,  the  matrix  could  be  outside  of  the  constraint 
region.  It  is  not  clear  from  this  representation  whether  a  covariance  matrix  is  in  the  interior  or 


9A  Vandermonde  vector  is  a  (complex)  N-ve ctor  such  that  d,  =  [a9  a,1  ••■af'-1]*.  In  this  case, 
ja,l  =  c.  Uniform  linear  arrays  have  the  property  that  the  direction  vectors  are  a  scalar  multiple 
of  Vandermonde  vectors. 
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exterior  of  the  constraint  region,  and  other  methods  would  be  needed  to  determine  whether  the 
matrix  was  realizable  for  a  particular  array  geometry.  This  inadequacy  of  this  representation  will 
be  illustrated  in  a  later  section. 

5.4  Inclusion  of  Receiver  Noise 

We  have  been  investigating  the  space  of  covariance  matrices  that  is  generated  by  a  spatial 
distribution  of  energy.  It  can  also  be  assumed  that  there  is  an  additional  additive  noise  term  that 
is  attributed  to  receiver  noise.  The  noise  covariance  matrix  attributed  to  this  term  will  add  to 
the  output  covariance  matrix  for  spatially  distributed  noise.  There  will  then  be  a  constraint  space 
for  the  output  covariance  matrix  that  is  the  cone  of  the  spatially  generated  noise  with  the  apex 
of  the  cone  shifted  awray  from  the  origin  by  the  covariance  matrix  of  the  receiver  noise.  If  the 
noise  covariance  matrix  is  unknown,  then  the  space  of  possible  output  covariance  matrices  will  be 
the  constraint  space  for  the  spatially  generated  noise  shifted  over  all  possible  shifts.  Unless  the 
noise  covariance  matrix  is  structured,  then  the  resulting  covariance  matrix  cannot  be  structured. 
Fortunately,  the  covariance  matrices  of  the  additive  receiver  noise  are  usually  structured. 

If  the  receiver  noise  is  assumed  to  be  independent  but  unknown,  the  dimension  of  the  con¬ 
straint  space  may  be  changed,  resulting  in  additional  terms  that  must  be  added  to  the  basis  func¬ 
tions  in  order  to  represent  the  constraint  space  or  a  covariance  matrix.  The  basis  functions  that  can 
be  added  to  account  for  this  are  the  outer  products  of  the  elementary  vectors.10  If  the  covariance 
matrix  due  to  the  spatially  distributed  noise  sources  can  be  represented  as  cqdjd,1,  the  entire 
covariance  matrix  can  be  represented  by 

L  Z.-.V 

R  =  tfidjdj  +  Y,  (125) 

t=l  i=L- t-1 

If  the  set  of  all  shifts  is  within  the  constraint  space  (e.g.,  due  to  independent  and  identically 
distributed  noise  sources  for  uniform  arrays),  then  the  combined  noise  covariance  matrix  will  retain 
the  structure  of  the  spatial  covariance  matrix,  and  the  representations  for  the  constraint  space  and 
covariance  matrix  will  apply  to  the  combined  constraint  space. 

Often  the  receivers  are  fully  characterized,  and  the  receiver  noise  variances  can  be  assumed 
known.  This  information  can  be  utilized  in  the  estimation  procedure  for  the  covariance  matrices. 

5.5  Discussions  and  Examples  for  Two  Sensors 

Examples  of  constraint  spaces  and  the  representation  of  some  elements  of  these  constraint 
spaces  will  now  be  shown.  The  examples  are  for  2  by  2  covariance  matrices.  The  arrays  that  are 


10An  elementary  vector  e,  is  zero  except  for  element  i,  which  is  unity. 
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used  to  form  the  matrices  are  assumed  to  be  narrowband  arrays,  and  the  array  element  responses 
are  identical.  The  covariance  matrices  that  result  for  these  arrays  are  Toeplitz  because  the  diagonal 
elements  will  be  equal,  and  thus,  the  dimension  of  the  constraint  space  is  3. 

Figures  will  be  shown  to  illustrate  the  constraint  spaces.  These  figures  are  generated  by 
noting  the  covariance  matrices  that  result  from  impulsive  spatial  spectra  or  a  discrete  interference 
source.  A  single  impulse  in  the  spatial  spectra  will  give  rise  to  a  covariance  matrix  that  is  on  the 
boundary  of  the  constraint  space.  This  matrix  is  a  weighted  outer  product  of  a  direction  vector 
with  itself.  In  this  case,  the  direction  vector  is  simply  a  vector  of  phase  shifts;  i.e.,  d  =  [le^']7". 
The  covariance  matrices  that  are  generated  by  a  spectral  impulse  are  matrices  of  the  form 


fin  fine  je' 
i?neJ0'  fin 


(126) 


From  Chapter  2,  for  a  narrowband  source  at  a  spatial  direction  6S ,  weight  a.  and  wavelength  A  and 
with  the  distance  between  sensors  s,  then 

R12  =  —  eJ2«cos«M/A  =  R  jo,  (127) 

‘  2tt  ' 

Varying  the  direction  of  arrival  of  the  source  9S  will  vary  9e  over  a  range  that  is  dependent  upon  the 
source  wavelength,  the  spacing  of  the  array  elements  and  the  possible  variations  in  the  direction  of 
arrival. 

The  constraint  space  is  then  generated  as  a  convex  combination  of  these  boundary  elements. 
The  covariance  matrices  resulting  from  a  single  impulse  in  the  spatial  spectra  are  singular  since  the 
determinate  of  the  matrix  is  0. 

Two  types  of  figures  are  used  to  illustrate  the  constraint  space  and  the  representations. 
The  first  figure  is  a  perspective  view  of  an  open  cone,  which  is  used  to  illustrate  the  covariance 
matrices  that  are  possible  for  the  constraint  spaces.  This  figure  is  plotted  in  the  coordinate  system 
[.ftii  5?(fii2)  3(fii2)jT  The  second  figure  is  a  cross  section  of  the  cone  for  /?n  =  R22  =  1-  The 
coordinate  axis  has  been  rotated  by  90  degrees  for  clarity.  The  reference  9e  =  0  is  along  the  3?(i?i2 ) 
axis. 


5.5.1  Covariance  Matrices  Possible 

The  first  example  will  be  for  a  narrowband  receiver  array  wrhere  the  two  elements  are  spaced 
one-half  of  a  wavelength  apart.  Both  of  the  elements  are  assumed  to  have  the  same  spatial  response 
and  gain.  For  this  array,  the  set  of  possible  covariance  matrices  is  an  open  cone,  as  shown  in  Figure  9. 
This  is  a  direct  result  of  the  integral  that  defines  the  space  of  covariance  matrices.  Substituting 
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the  element  separation  into  Equation  (127),  the  exponent  is  jn  cos(6a).  Varying  9,  over  the  domain 
0  <  9S  <  7r  results  in  the  range  for  the  electrical  angle  -tr  <  9e  <  n. 


K{«12} 


3{fl12} 


flll- ^22  S{fil2} 


Figure  9.  Constraint  space  for  one-half  wavelength  spacing. 


The  second  example  will  be  for  a  narrowband  array  where  the  array  elements  are  placed 
three-eighths  of  a  wavelength  apart.  The  set  of  possible  covariance  matrices  for  this  array  is  also  an 
open  cone  but  has  a  hyperplane  for  a  portion  of  the  boundary  rather  than  having  the  continuous 
support  of  the  previous  example.  This  is  shown  in  Figure  10.  This  is  a  result  of  the  mapping  from 
the  physical  angle  of  the  source  to  the  electrical  angle.  Proceeding  as  for  the  previous  example,  the 
exponent  for  the  R 12  term  is  j|7rcos(0j).  Varying  9S  over  the  domain  0  <  9S  <  7r  results  in  the 
range  for  the  electrical  angle  -\ir  <9e<  |tt. 

The  range  of  9e  here  has  been  restricted  by  the  spacing  of  the  array  elements.  This  restriction 
on  9e  could  also  be  achieved  if  the  domain  of  the  angles  of  incidence  of  possible  interferes  is 
restricted. 


5.5.2  Representations 

The  first  representation  is  to  approximate  the  constraint  set  bv  the  use  of  the  finite  approx¬ 
imation  to  the  covariance  integral.  The  cone  formed  by  the  finite  basis  is  shown  in  Figure  11  for 
four  terms  in  the  sum  as  well  as  for  eight  terms.  The  terms  in  the  summation  are  spaced  such  that 
the  covariance  lags  resulting  from  this  spacing  are  uniformly  spaced  in  electrical  angle. 
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*{Rn} 


R11.R22  ^{^12} 


Figure  10.  Constraint  space  for  three-eighth  wavelength  spacing. 


It  is  assumed  that  a  particular  covariance  matrix  is  within  the  constraint  set  formed  by  the 
finite  approximation.  The  number  of  terms  in  the  summation  used  to  generate  this  matrix  can  be 
reduced  using  the  reduction  theorem.  For  this  example,  there  are  several  polyhedral  cones  that 
can  be  formed  containing  the  covariance  matrix.  This  is  shown  in  Figure  12  for  two  of  the  cones 
that  could  result.  The  intersections  of  the  polyhedral  cone  with  the  boundary  of  the  constraint  set 
indicate  the  frequencies  that  are  used  to  represent  the  covariance  matrix.  This  example  illustrates 
the  non-uniqueness  of  the  spatial  spectrum.  Both  of  these  representations  of  the  covariance  matrix 
are  maximal  representations. 

Figure  13  illustrates  a  hyperplane  that  intersects  the  constraint  space  and  has  the  indicated 
covariance  matrix  as  an  element.  There  are  many  hyperplanes  that  have  this  property.  The 
hyperplanc  shown  is  the  hyperplane  that  is  indicated  by  the  Grenander  and  Szego  theorem  of 
Caratheodorv.  This  represents  the  covariance  matrix  as  a  weighted  sum  of  the  vectors  labeled  .4 
and  B  in  the  three-dimensional  projection  of  Figures  13  and  14  or  as  a  convex  combination  of  the 
points  labeled  .4  and  B  in  the  cross  section.  For  the  array  with  three-eighths  wavelength  spacing, 
the  vector  (point)  B  is  not  a  member  of  the  constraint  set.  and  there  is  not  a  physical  angle  for 
a  source  that  would  generate  this  vector.  The  spectral  point  indicated  is  not  valid  for  this  array 
geometry. 
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X{Rn}  R{R  12} 
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Figure  11.  Finite  sum  approximation  and  resulting  constraint  space 
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Figure  13.  Hyperplane  representation. 
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Figure  14-  Invalid  Caratheodory  representation. 
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5.6  Conclusion 


Representations  for  the  covariance  matrices  using  a  finite  number  or  parameters  have  been 
introduced  in  this  chapter.  These  representations  and  the  properties  of  the  constraint  set  will  be 
used  to  form  estimates  of  covariance  matrices  for  use  in  adaptive  signal-processing  algorithms. 
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6.  FIXED-BASIS  ESTIMATOR 


6.1  Introduction 

The  first  approach  that  will  be  used  to  form  joint  estimates  of  the  structured  mean  and 
covariance  matrix  is  to  approximate  the  integral  expression  that  characterizes  the  constraint  space 
using  the  finite  sum  representation  discussed  in  the  previous  chapter.  We  will  then  derive  an 
estimator  for  the  mean  and  covariance  matrix  based  on  this  spectral  representation  of  the  constraint 
space.  The  resulting  estimator  can  be  restricted  so  that  only  positive-definite  estimates  of  the 
covariance  will  result.  The  finite  sum  approximation  may  introduce  a  bias  for  the  covariance 
estimates  since  not  all  spatially  constrained  covariance  matrices  can  be  represented  by  the  finite 
sum. 

For  simplicity  the  system  is  assumed  to  be  narrowband.  The  covariance  matrix  will  have  the 

form 

f  S{aj0, 9.  o)d(u;o.  0)d\jjo,&, 0)  sin  Qdddd)  diag(r/j  ...  rj\)  (128) 

2?r  Je.o 

The  dependence  on  will  be  suppressed.  Here,  the  terms  that  correspond  to  receiver  noise 
are  assumed  unknown.  Later,  the  modifications  necessary  when  these  terms  are  either  known  or 
assumed  equal  to  each  other  will  be  discussed. 

The  data  vectors  are  ordered  such  that  the  first  G  sample  vectors  are  assumed  to  be  non¬ 
zero  mean,  and  the  next  K  sample  vectors  are  assumed  to  be  zero  mean.  The  sample  vectors  are 
assumed  to  be  independent  and  are  arranged  in  a  data  matrix  Z  =  [zi  •  •  •  *.p\ ,  P  —  G  +  K.  The 
mean  is  of  the  form  bkd(6,o).  corresponding  to  an  unknown  amplitude  6*  multiplying  a  known 
direction  vector  parameterized  by  the  direction  of  arrival  9.0-  The  multivariate  density  functions 
underlying  all  of  these  samples  are  assumed  to  share  the  same  covariance  matrix.  The  log-likelihood 
for  the  data  is  then  given  by 

l  =  - NP log  7r  -  Flog  |jR|| 

a  p 

-Y,(zk~bkd(9.o))1R-\zk-bkd(lo))-  £  (129) 

We  would  like  to  find  a  maximum-likelihood  estimator  for  the  mean  and  the  covariance  under  this 
model. 

In  this  section  the  covariance  estimates  will  be  found  using  the  M- term  finite  sum  approxi¬ 
mation  for  the  covariance  integral  in  Section  5.3.2. 
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(130) 


M 

R  =  5Z  5(^m,om)d(0m,<?rn)dt(em,<?m)  +  d\ag(r)l . . .  T)%)  =  DSD1 

m=l 


D  is  an  Nx(N  +  M)  matrix  formed  from  an  NxM  array  of  steering  vectors  with  the  identity  matrix 
appended. 


D  = 


d(#i,  d>i)  ■  ■  ■  d(8\f,  0a/):I.v 


(131) 


and  S  is  a  diagonal  matrix  with  non-negative  diagonal  entries. 

6.2  Derivation  of  the  EM  Algorithm 

A  closed-form  expression  for  the  maximum-likelihood  estimator  for  the  covariance  has  not 
been  found;  however,  the  Expectation-Maximization  (EM)  algorithm  [47]  can  be  used  to  derive  an 
iterative  algorithm  where  each  iteration  will  remain  within  the  constraint  region  and  each  iteration 
will  not  decrease  in  likelihood. 

The  EM  algorithm  consists  of  the  following  two  steps. 

1.  E-step:  Given  a  complete  data  set  Y.  that  if  known  would  uniquely  determine  the 
observed  data,  compute 

E{/cd(Z;0)!Y,©m}  (132) 

lcd\ Z:  0)  here  is  the  complete  data  log-likelihood,  not  the  log-likelihood  for  the  mea¬ 
sured  data. 

2.  M-step:  Find 

0m*1  =  argmax  E{lcd(Z:  ©)|Y;  0m}  (133) 

© 

The  EM  algorithm  on  be  shown  to  produce  a  sequence  of  estimates  for  which  the  likelihood  is 
non-decreasing. 

For  the  EM  algorithm  two  data  spaces  must  be  defined.  The  complete  data  space  is  a 
hypothetical  data  space,  and  we  will  choose  M  4-  fV-length  vectors  y  arranged  in  a  matrixY  = 
[yi  •  ■  •  y p]  such  that 

Z  =  DY  (134) 

The  covariance  of  the  columns  of  Y  will  be  a  diagonal  matrix  E  with  the  first  M  diagonal  terms 
being  the  discretized  spectrum  S  and  the  next  N  diagonal  terms  corresponding  to  receiver  noise. 
The  incomplete  data  space  is  our  observed  data  Z.  The  matrix  D  will  be  ordered  so  that  the 
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direction  vector  corresponding  to  the  mean  is  the  first  column  of  the  matrix  D  and  will  be  denoted 
dj.  The  vectors  y*  that  correspond  to  the  non-zero  mean  zk  will  consist  of  a  noise  term  plus  an 
additive  signal  term  that  is  a  multiple  of  the  first  elementary  vector 


yk  =  bkex  +  m 


(135) 


with  dj  =  Dei. 

The  complete  data  log-likelihood  is 


led  =  -(A/  +  N)P  log  n-P  log  ||£|| 

G  P 

-^(yg-bket)^~l(yk-bkei)-  ]T  yl-E^y* 

fc  =  l  k=G+ 1 


(136) 


If  the  diagonal  entries  of  £  are  represented  as  a\  .  . .  the  log-likelihood  can  then  be  written 

m  -  .v 


cd 


(M  +  .V)Plog7T  -  P  ^  log<7_,  - 


J=1 
M-.V  G 


o  ,  ,1  m-.n  o  |,r  1  Af-.V  P 

El  *  lit  —  Ofc  ^  ^  _  W'  >  jk 

d\  Z_r  or  Z-r  (7 

fc  =  l  1  j  =  2  fc=l  J  j  =  l  *  =  G-1  J 


(137) 


The  definition  of  Y  has  been  used  to  avoid  double  subscripts  since  Yjk  =  (y k)j. 

Define  B  =  j&i  •  •  •  6cj.  The  expected  value  of  the  log-likelihood  can  be  found  (the  E-step)  to 


be 


E{lcd'±p.Bp.Z}  = 


.u-.v 


( M  .VlPlogTr  -  P  logoy  -  ^ 


E{|l'u  -  6*i2  !£p,  Bp.  Z} 


j  =  i 


*=i 


<y  i 


A^V  £  E{!y>j2i£p.Bp.Z}  _  A^V  £  E{lljt|2  ]£p.  Bp.  Z} 

j  =  2  k=  1  j  =  l  k=G+l  J 

The  third  term  in  Equation  (138)  can  be  expanded  as  follows: 

E{IVu  -  bk\2  i£p.  Bp.  Z}  =  E{iVu|2  ]£p.  Bp.  Z}  -  2Re(6;E{Yifc  |£p.  Bp.  Z})  +  !6tj2 

=  E{T'ui2!£p.Bp.Z} 

-iE{V'i* [£p.  Bp.  Z}|2  +  1E{>T*  |EP.  Bp,  Z}  -  bk\2 


(138) 


(139) 
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The  conditional  expectations  can  be  written  as 


E{y*|Ep,Bp,Z}  =  ^1 +  £^11?  ‘(z*-6£di)  (140) 

and 

E{|}>!2;Sp,Bp  Z}  =  [Sp-EpDtRp'1DEp]JJ  +  jE{(rjfc|Ep  BP,Z}|2  (141) 


for  the  data  that  have  non-zero  means,  and 


E{j>'jfc|2|£p,Bp.Z}  =  [EP-S^RP  ‘dEp  +  SpDfRp  ^z^Rp  (142) 


for  the  data  that  are  zero  mean. 

Next,  the  arguments  that  maximize  Equation  (138)  (the  M-step)  can  be  found.  Maximizing 
this  first  with  respect  to  the  unknown  b k  terms  yields 


&T1  =E{rlfc|Sp.Bp,Z}  =  bp  +  [EpDtRp  X(z*-6gdi)] 


(143) 


Substituting  this  in,  we  next  find  the  a3  that  maximize  Equation  (138). 


*r‘  =  < 


p 


k= i 
K 


P 

k=K  i-l 


^E{,}-fc|2|Ep,Bp  Z}+  Y.  E{|}jfci2!Ep, Bp,Z} 


J  >  1 


£E{!ylfc|2|Sp.B.Z} 


Lfc=l 


Y  E{iyu!2,Sp.Bp.Z}- ;E{l',;Sp.Bp.Z}| 

k=K~  1 


The  resulting  iterations  are 

±pj'  -  [sp  -  SpDf  (Rp-1SpRp_1  -  Rp_1)  D£p 


J  =  1 


RP~‘  =  DEp+1Dt 


(144) 


(145) 

(146) 


and 
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fcT1  =%+  [spDf  RP 


(147) 


where  the  sample  covariance  matrix  is 


5Z  z*zl 

'  G  P 

J2(zk  -  ^di)(zfc  -  ^di)f  +  ^  Zfczj. 


J  =  1 

j  # 1 


(148) 


If  the  true  covariance  were  known,  the  EM  iterations  would  not  be  needed  since  closed-form 
expressions  for  the  unique  maximizing  b k  are  well-known.  This  observation  leads  us  to  modify  the 
EM  iterations  so  that  the  likelihood  with  respect  to  the  bk  is  directly  maximized  in  each  iteration. 
The  expressions  for  the  bk  that  will  accomplish  this  are 


if  =  4  *-> 

dj  RP  di 


(149) 


Modifying  the  iteration  rule  for  the  bk  allows  us  to  combine  the  iteration  rules  for  all  of  the 
elements  of  £  into  a  single  iteration  rule  using  the  sample  covariance  matrix  in  Equation  (148)  for 
j  ^  1.  The  non-zero  mean  data  will  add  to  the  sample  covariance  such  that  the  iterations  on  £  are 
unaffected  by  combining  the  iteration  rules.  This  can  be  shown  bv  evaluating  products  that  would 
then  appear  in  the  iteration  rule  for  j  =  1. 

dt1RP'1(zfc  -  ftgd,)  =  d\W>-\k  -  djRP^dx  dj  ^  _*k 

dj  Rp  dj 

=  0  (150) 

It  can  then  be  seen  that  the  iterations  are  unaffected  by  combining  the  rules.  This  also  shows  that 
the  mean  terms  would  not  change  for  the  next  iteration  of  the  formal  EM  algorithm.  This  term 
enters  directly  into  the  update  equation  for  the  mean.  By  making  this  modification  to  the  formal 
EM  algorithm,  the  direction  vector  corresponding  to  the  mean  need  not  be  one  of  the  columns  of 

D 

The  algorithm  that  results  from  this  modification  is  a  generalization  of  the  EM  algorithm.  It 
has  the  property  that  the  likelihood  produced  by  the  sequence  of  estimates  is  non-decreasing,  and 
for  each  iteration  the  increase  in  likelihood  is  lower  bounded  by  the  increase  in  likelihood  of  the 
formal  EM  iteration.  This  can  be  seen  by  examining  the  estimate  of  the  mean  produced  by  the 
formal  EM  iterations.  If  the  likelihood  is  maximized  with  respect  to  the  mean  conditioned  on  the 
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current  estimate  of  R.  this  would  also  be  the  estimate  for  the  modified  algorithm.  If  the  likelihood 
is  not  maximized  fully  with  respect  to  the  mean,  replacing  the  formal  EM  mean  estimate  with  the 
modification  will  result  in  an  increase  in  likelihood. 

The  generalized  EM  iterations  are 

±pj1  =  jsp  +  SPD1  (RP_1SpRp_1  -  Rp_1)  DEP]^  ,  (151) 


Rp+1  =  DEp+1Df 


(152) 


and 


d|  R  zk 
d{  R  ‘d! 


where  the  sample  covariance  matrix  is 


G 

^ di)(2fc  -  fydrf 

k  =  l 


E  z*zI 


(153) 


(154) 


If  all  of  the  samples  are  assumed  to  be  zero  mean,  then  this  EM  algorithm  is  similar  to  the 
algorithm  of  Miller  and  Snyder  [26]  and  Moulin  et  al.  [48]  where  it  was  proposed  for  Toeplitz 
constrained  covariances. 

The  iterations  for  the  spectral  estimates  take  the  form  of  a  gradient  ascent  algorithm  since 
the  diagonal  terms  of 

Df  ^Rp_'spRp_1  -  Rp_1j  D  (155) 

are  the  gradients  of  the  log-likelihood  with  respect  to  each  of  the  elements  of  E.  At  stationary 
points  of  the  iterative  algorithms,  the  gradient  will  either  be  zero  for  the  elements  of  E  that  are 
non-zero  or  will  be  negative  or  zero  for  elements  of  E  that  are  zero.  The  estimates  of  the  elements 
of  E  will  never  be  negative  as  they  are  found  in  each  iteration  as  conditional  correlations.  For 
these  reasons  the  stationary  points  of  the  iterative  algorithm  coincide  with  local  maxima  of  the 
likelihood. 

If  the  receiver  noise  variances  are  assumed  to  be  equal,  then  the  algorithm  is  modified  by 
replacing  each  of  these  terms  by  their  average  in  each  iteration.  This  is  a  straightforward  modifi¬ 
cation  that  appears  in  the  maximization  step  of  the  derivation.  If  the  receiver  noise  variances  (or 
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similarly  the  power  level  of  the  interference  from  any  direction)  are  known,  then  these  terms  do 
not  need  to  be  estimated  and  can  be  replaced  in  each  iteration  bv  the  known  values. 

6.3  Existence  of  Positive-Definite  Solutions 

In  order  to  determine  if  a  positive-definite  solution  exists  to  the  problem  characterized  above, 
we  will  first  assume  that  the  receiver  noise  terms  are  known  and  equal  to  zero.  Additionally,  it  will 
be  assumed  that  the  covariance  matrix  is  estimated  from  data  that  are  known  to  be  zero  mean. 

The  conditions  necessary  for  the  likelihood  to  be  unbounded  above  were  stated  for  a  class  of 
covariance  matrices  R  satisfying  certain  restrictions  [49].11  Under  these  restrictions,  the  likelihood 
will  be  unbounded  above  if  for  a  set  of  independent  observations  Z  —  {z\  -  ■  ■  z x},  a  singular  matrix 
Ro  exists  such  that  zjt  is  in  the  range  space  of  Ro  for  k  —  1  •  •  •  K. 

Our  problem  is  parameterized  such  that  the  restrictions  on  the  class  of  covariance  matrices 
R  are  satisfied.  The  likelihood  will  be  bounded  above  with  probability  1  since  the  requirement 
that  a  singular  matrix  exists  with  the  data  in  the  range  space  is  satisfied  with  probability  0,  even 
for  A"  =  1.  There  is  a  finite  number  of  reduced  dimensional  subspaces  corresponding  to  singular 
covariance  matrices,  and  the  probability  that  any  z *  lies  in  one  of  these  subspaces  is  0. 

Our  attention  will  now  turn  to  the  case  where  all  of  the  data  are  assumed  to  contain  means 
that  must  be  estimated.  A  maximum-likelihood  estimate  for  this  data  set  may  not  exist.  If  there 
exists  a  maximal  proper  subspace  whose  V  -  I  basis  vectors  in  D  are  orthogonal  to  dj,  then  the 
likelihood  is  unbounded  above,  and  the  maximum  likelihood  estimate  of  R  does  not  exist.  The 
existence  of  this  subspace  is  a  deterministic  problem,  and  it  sets  conditions  on  the  vectors  that  are 
the  columns  of  D. 

Non-existence  of  the  maximum-likelihood  estimate  can  be  shown  by  examining  the  likelihood 
as  certain  of  the  parameters  of  the  likelihood  are  allowed  to  become  arbitrarily  small.  Let  the  matrix 
C  be  the  matrix  formed  from  the  weighted  sum  of  outer  products  of  the  orthogonal  subspace  basis, 
di  will  be  an  eigenvector  of  C  with  associated  eigenvalue  0.  Let  the  other  eigenvalues  of  C  be 
denoted  by  62  ■  ■ .  <5.v  >  0.  Now  form  a  covariance  matrix 

R  =  C^Mid{  (156) 


with  inverse 


R  1  =  C  *  -  y-didt 

Cl 


(157) 


nThis  paper  was  written  with  application  to  real  symmetric  matrices.  It  is  not  difficult  to  show 
that  the  conditions  are  true  for  complex  matrices  as  well. 
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C  *  is  the  Moore- Penrose  pseudoinverse  of  C  [36],  di  is  an  eigenvector  of  R  and  R  l. 

Now  substitute  this  covariance  matrix  into  the  log-likelihood  functional  for  a  single  sample 
and  investigate  the  effect  of  b\  becoming  arbitrarily  small.  For  a  single  sample,  the  log-likelihood 
with  constant  terms  removed  is 


l  =  -  log  |jR||  -  (z  -  6di)tR  1(z  ~  bdi) 


(158) 


Substituting  the  maximum-likelihood  estimate  b  = 
result  is 


and  the  matrix  R  and  simplifying,  the 


v 

Z  =  -£ log<5, -zfCr*z  (159) 

i  =  l 

As  approaches  zero,  the  log-likelihood  l  approaches  infinity;  thus,  the  likelihood  is  unbounded 
above  and  a  maximum- likelihood  estimate  will  not  exist.  This  will  hold  true  for  an  arbitrarily  large 
number  of  data  samples  provided  a  mean  is  to  be  estimated  for  each  sample. 

In  the  discretization  of  the  integral  equation  for  the  covariance  matrix,  it  was  assumed  that 
the  discretization  would  adequately  represent  the  resulting  covariance  matrix.  In  order  to  form 
this  representation,  a  large  number  of  spectral  points  and  direction  vectors  may  be  needed.  The 
existence  of  the  condition  allowing  the  likelihood  to  be  unbounded  can  be  determined  by  examining 
whether  .V  - 1  columns  of  D  are  orthogonal  to  di .  Restricting  the  basis  vectors  in  this  manner  is  not 
appealing:  rather,  the  data  model  should  be  refined.  It  is  known  that  the  receiver  noise  variances 
are  non-zero  and  can  often  be  estimated  accurately.  Including  the  receiver  noise  variance  in  the 
covariance  model  will  lower  bound  the  minimum  eigenvalue  of  the  covariance  matrix  preventing  a 
singular  covariance  estimate.  The  likelihood  will  then  be  bounded  above. 

6.4  Asymptotic  Properties 

The  asymptotic  behavior  of  the  estimates  will  indicate  an  upper  bound  on  the  performance 
of  any  beamformer  or  detector  that  uses  the  estimates.  We  are  interested  in  the  behavior  for 
two  asymptotes:  for  the  amount  of  data  increasing  and  for  the  number  of  samples  in  the  discrete 
approximation  increasing. 

First,  the  behavior  of  the  covariance  estimates  when  the  number  of  zero  mean  data  approaches 
infinity  will  be  examined.  The  sample  covariance  matrix  will  converge  to  the  true  covariance  matrix 
R<r«<?  £  fZj  with  probability  1.  If  this  matrix  can  be  represented  by  the  discretized  model,  then 
this  will  be  the  maximum-likelihood  estimate,  and  this  is  a  stable  point  of  the  EM  algorithm.  The 
estimator  is  then  asymptotically  unbiased. 
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If  the  true  covariance  is  near  singular,  then  there  may  not  be  a  representation  of  the  true 
covariance  in  terms  of  the  discretized  spectrum,  and  an  asymptotically  biased  estimator  results. 
All  implementations  of  this  estimator  must  be  formed  from  a  finite  number  of  spectrum  samples 
and  will  not  represent  all  possible  covariance  matrices.  The  asymptotic  bias  of  the  estimator  will 
result  in  a  loss  in  likelihood  compared  to  the  maximum,  impacting  detectors  that  are  based  on 
likelihood.  It  would  also  be  expected  that  the  loss  in  signal-to-noise  ratio  for  beamformers  that  use 
the  covariance  estimate  would  approach  zero  as  the  amount  of  data  becomes  large.  Because  the 
true  covariance  cannot  be  represented,  the  optimal  weight  vector  will  not,  in  general,  be  calculated, 
and  a  loss  in  signal-to-noise  ratio  will  occur. 

6.5  Brief  Simulation  Results 

The  performance  gain  of  an  adaptive  beamformer  based  on  this  constrained  estimator  will 
be  used  to  show  that  adaptive  signal  processing  using  constrained  covariance  estimates  is  viable. 
The  loss  in  signal-to-noise  ratio  for  adaptive  beamformers  using  constrained  and  unconstrained 
estimates  of  the  covariance  matrix  will  be  compared.  This  loss  factor 


[dtR._1di2 

d+R-1ddtR-IRR'1d 


(160) 


is  a  random  variable,  and  the  histogram  of  the  loss  factor  for  constrained  covariance  estimates  will 
be  compared  with  the  analytic  expression  for  the  density  function  when  unconstrained  estimates 
are  used.  The  histograms  are  found  using  Monte  Carlo  simulation. 

First,  the  performance  is  shown  when  the  true  covariance  is  within  the  polyhedral  cone  of  the 
finite  basis  functions.  The  array  simulated  here  is  a  4-element  minimum  redundancy  array  with 
A  =  .V,  This  array  and  the  simulations  will  be  discussed  further  in  Chapter  9.  Figure  15  illustrates 
the  performance  that  can  be  obtained  using  structured  estimates  compared  to  the  unconstrained 
covariance  estimate.  The  line  labeled  ML  is  the  density  function  of  the  signal-to-noise  ratio  loss 
factor  when  the  unconstrained  estimate  of  the  covariance  is  used  in  an  adaptive  beamformer.  The 
line  labeled  CML  is  an  estimate  of  the  loss  factor  density  function  when  an  adaptive  beamformer 
is  based  on  the  constrained  covariance  estimate.  This  is  the  histogram  for  200  realizations.  The 
increase  in  the  mean  of  the  loss  factor  indicates  that  a  dramatic  increase  in  the  performance  of 
adaptive  signal  processors  could  be  expected  when  using  the  constrained  estimates. 

In  Figure  16  the  effects  of  an  insufficient  number  of  terms  in  the  finite  sum  are  shown.  The 
true  covariance  here  is  outside  of  the  assumed  polyhedral  cone.  The  density  function  shows  a  high 
loss  in  signal-to-noise  ratio  when  compared  to  the  constrained  covariance  estimate  with  an  adequate 
number  of  terms  in  the  finite  sum.  A  beamformer  based  on  this  covariance  estimate  would  also 
show  a  loss  when  compared  to  the  beamformer  using  an  unconstrained  covariance  estimate. 

In  the  next  chapter,  an  alternative  parameterization  of  the  covariance  will  be  introduced. 
This  parameterization  is  based  on  the  representation  theorems  introduced  in  Chapter  5  and  will  not 
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Figure  15. 


Figure  16. 
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Signal-to-noise  ratio  loss  factor  density.  Adequate  number  of  terms. 
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require  the  large  number  of  terms  that  may  be  required  to  approximate  adequately  the  covariance 
integral. 
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7.  VARIABLE  BASIS  ESTIMATOR 


The  reduction  theorem  and  the  Caratheodory  representation  theorem  discussed  in  Chapter  5 
state  that  any  member  of  an  L-dimensional  convex  set  can  be  represented  as  the  convex  combi¬ 
nation  of  L  +  1  elements  of  that  set.  The  large  number  of  spectral  weights  used  in  the  previous 
chapter  to  induce  the  constraint  is  mathematically  unnecessary  and,  because  of  the  large  number 
of  calculations  needed  to  calculate  the  weights,  is  undesirable  from  an  engineering  perspective.  In 
this  chapter,  a  method  is  proposed  that  can  be  used  to  form  joint  estimates  of  the  structured  means 
and  covariance  matrices  by  using  the  maximal  representation.  The  proposed  method  is  capable  of 
representing  all  covariance  matrices  satisfying  the  spatial  constraint.  The  fixed-basis  finite  repre¬ 
sentation  used  in  Chapter  6  is  unable  to  accomplish  this.  The  reduction  in  the  number  of  spectral 
terms  that  must  be  estimated  reduces  the  per-iteration  computational  demand  and  increases  the 
rate  of  convergence  of  the  resulting  covariance  estimates. 

The  Caratheodory  representation  theorem  can  be  viewed  as  stating  that  there  is  an  inscribed 
polyhedral  cone  formed  from  L  ~  1  bases  containing  an  estimated  covariance  matrix.  This  cone 
cannot  be  formed  a  priori  because  the  covariance  matrix  may  take  on  values  anywhere  within  the 
constraint  region.  For  each  estimated  covariance  matrix,  a  polyhedral  cone  must  be  found  that 
contains  the  estimate  in  the  interior.  To  determine  this  cone,  the  direction  vectors  used  to  form 
the  basis  functions  must  be  allowed  to  change  by  varying  the  directions  of  arrival.  This  is  very 
similar  to  the  stochastic  direction  of  arrival  problem  where  it  is  assumed  that  the  spatial  source 
contribution  to  the  covariance  estimate  lies  on  a  boundary  of  the  constraint  space.  The  angles 
parameterizing  the  basis  functions  for  this  hyperplane  are  estimated  to  determine  the  estimated 
directions  of  arrival. 

It  would  be  desirable  to  use  current  direction-finding  algorithms  if  they  could  be  extended  to 
fit  our  requirements.  Unfortunately,  these  methods  suffer  from  shortcomings  that  limit  their  useful¬ 
ness  for  our  purposes.  Eigen-structure  or  subspace  methods,  such  as  multiple  signal  classification 
(MUSIC),  require  that  the  number  of  interference  sources  be  limited  to  A  -  1  sources  and  full  rank 
sample  covariance  matrices.  Previous  stochastic  direction-finding  algorithms  based  on  maximizing 
likelihood  are  known  to  perform  poorly  when  there  are  many  interference  sources  and  the  sources 
are  spatially  closer  together  than  one-quarter  of  a  beamwidth  [50  .  Additionally,  the  stochastic 
direction-finding  algorithms  require  a  computationally  intensive  search  over  the  possible  directions 
of  arrival  of  the  interference  sources  for  each  iteration. 

A  new  procedure  based  on  the  estimator  of  Chapter  6  but  which  does  not  fix  the  parameters 
of  the  liases  will  be  introduced.  This  estimation  procedure  allows  the  estimation  of  any  covariance 
matrix  that  is  within  the  constraint  space.  The  likelihood  is  non-decreasing  for  the  sequence 
of  estimates  provided  by  the  algorithm,  and  stable  points  of  the  resulting  algorithm  satisfy  the 
necessarv  conditions  required  for  a  maximum-likelihood  -stimate. 


7.1  Preliminaries 


In  this  section  we  will  briefly  digress  in  order  to  understand  some  items  that  will  be  used  in 
the  proposed  estimation  algorithms.  Three  topics  of  interest  along  with  the  Caratheodory  theorem 
are  used  to  justify  the  proposed  algorithm.  The  first  topic  is  how  the  likelihood  surface  will  affect 
the  constrained  covariance  estimates.  The  second  topic  is  a  discussion  of  how  the  discretization 
used  in  Chapter  6  affects  convergence  of  the  estimate  for  certain  estimated  covariance  matrices.  A 
maximum-likelihood  spectral  estimator  is  then  discussed. 

7.1.1  Likelihood  Surface 

The  likelihood  is  unimodal  over  all  covariance  matrices  for  a  non-singular  sample  covariance 
matrix  and  has  only  a  single  local  maxima.  If  the  estimate  moves  in  any  direction  away  from 
the  unconstrained  ML  estimate,  the  likelihood  is  a  monotonically  decreasing  function.  Forming 
covariance  estimates  Ra  =  StC  and  R;,  =  S  4-  aC  with  RQ,R{,  >  0  and  with  C  =  Ct  then,  for 
q  >  1,  /(Ra)  >  /( Rb)  with  equality  only  if  C  =  0.  S  is  the  sample  covariance  matrix  that  provides 
the  maximum  in  the  unconstrained  case. 

The  unimodality  of  the  likelihood  is  shown  in  the  following  proof.  S  and  C  are  similar  to 
diagonal  matrices:  therefore,  there  exists  a  transformation  A  that  will  simultaneously  diagonalize 
S  and  C  [281.  Define  the  diagonal  matrices  S  =  A-1SA  and  C  =  A~*CA  with  diagonal  elements 
S!  ...  i  v  and  cj  . . .  c\ .  This  transformation  is  linear  and  will  simultaneously  diagonalize  Ra  and 
Rf,.  Using  these  definitions  and  neglecting  constant  terms,  the  log-likelihoods  are  proportional  to 


l(Ra)  = 

i=l 


log(s,  t-  cn)  + 


S, 

s,  H-  C, 
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and 
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The  difference  in  log-likelihood  between  these  estimates  is 


/(Ra)  -/( Rfr)  = 


i=  1 


,  .  Si  +d  s,  s, 

log(  r - - )  -1-  V - r-  -  r - -  I 

St  +  ac,  s,  ~  ct  .s,  -  or,  j 


Applying  the  inequality  log  j-  <  J"  -  1.  then 
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with  equality  if  and  only  if  ac,  =  c,.  Simplifying  this  expression  yields 


Z(Ra)  - /(R6)  >  ]T  —  C‘(Q  ^ 


(S,  -+■  Ci)(3,  +  ac,) 


(165) 


All  of  the  terms  that  enter  into  this  sum  are  non-negative  for  a  >  1,  showing  the  unimodality 
of  the  likelihood  surface.  For  a  >  1,  the  equality  in  Equation  (165)  holds  only  if  all  c,  are  zero 
corresponding  to  C  =  0.  A  unimodal  surface  such  as  this  will  have  only  a  single  extremum. 

Unimodality  of  the  likelihood  surface  for  unconstrained  maximum-likelihood  estimates  does 
not  indicate  the  shape  of  the  likelihood  surface  over  more  restricted  covariance  spaces.  There  may 
be  ridges  in  the  unconstrained  likelihood  surface  even  though  it  is  unimodal.  If  a  constraint  set 
passes  through  several  of  these  ridges,  then  local  maxima  will  exist  and  the  likelihood  surface  will 
be  poorly  behaved. 

The  likelihood  is  continuous  for  positive-definite  covariance  matrices.  The  determinant,  ma¬ 
trix  inverse,  and  the  trace  functions  are  all  continuous  for  positive-definite  arguments. 

7.1.2  Constrained  Covariance  Estimates  and  Finite  Bases 

If  the  sample  covariance  is  within  the  constraint  space  (such  as  occurs  for  A'  —  oc).  then 
the  maximum-likelihood  covariance  estimate  will  be  the  sample  co%ariance  matrix.  Any  movement 
within  the  constraint  space  away  from  the  sample  covariance  will  result  in  a  corresponding  loss  in 
likelihood.  If  the  sample  covariance  matrix  is  not  w,  hin  the  polyhedral  cone,  then  the  likelihood  will 
not  reach  the  maximum  for  estimates  within  the  constraint  space.  When  a  constrained  maximum- 
likelihood  estimate  is  outside  of  the  polyhedral  cone  formed  by  the  finite  basis,  by  continuity  of 
the  likelihood  and  unimodality  over  all  matrices,  the  finite  basis  constrained  estimate  will  lie  on 
the  boundary  of  the  polyhedral  cone.  This  can  be  proven  by  contradiction.  Assume  that  there  is 
a  maximizer  interior  to  the  polyhedral  cone.  By  the  argument  of  Section  7.1.1.  the  likelihood  will 
be  increasing  along  the  line  from  the  interior  point  to  the  sample  covariance  exterior  to  the  cone; 
therefore,  the  interior  point  cannot  be  a  maximizer.  This  boundary  is  a  hyperplane  of  dimension 
at  most  L  -  1 . 

If  an  estimate  lies  on  a  hyperplane  boundary,  then  only  those  basis  elements  that  are  in  the 
subspace  defined  by  the  hyperplane  may  have  non-zero  weights.  If  any  set  of  L—  1  bases  are  linearly 
independent,  then  at  most  Z.  -  1  spectral  weights  will  be  non-zero.  This  condition  can  be  easily 
tested;  however,  it  is  not  altogether  useful  for  our  purpose. 


The  EM  estimator  performs  poorly  when  the  spectral  weights  are  approaching  zero,  that 
is.  when  the  covariance  estimate  is  approaching  a  hyperplane  boundary.  Examining  the  update 
equation  for  a  spectral  weight 

<r*+1  =  c^  +  a^d)  (rp_1SPRp_1  -  RP_1)d;^  .  (166) 

it  is  seen  that  this  can  be  written  as 

=  crf(1  +a*jA)  .  (167) 

The  term  A  here  is  the  gradient  of  the  likelihood  with  respect  to  ay  The  fractional  change  in 
the  value  of  the  spectral  estimate  will  be  given  by  rrj’A.  As  the  spectral  estimate  decreases,  then 
the  estimate  will  change  by  a  smaller  and  smaller  fraction  of  its  value  fc.  each  iteration.  The 
value  of  A  will  be  decreasing  as  well  when  the  estimate  moves  closer  to  a  maximum-likelihood 
estimate.  For  these  reasons,  the  spectral  estimates  will  approach  zero  very  slowly.  For  finite  word 
implementations  roundoff  error  in  calculating  the  update  will  result  in  a  spectral  estimate  that  does 
not  change  when  cr^A  is  smaller  than  the  machine  precision. 

We  would  like  for  any  estimate  within  the  constraint  region  to  be  represented  within  some 
finite  approximation  to  the  constraint  space.  When  the  estimate  is  moving  outside  of  the  polyhedral 
cone,  as  evidenced  by  the  solution  on  or  approaching  the  hyperplane  boundary,  then  the  cone  can 
be  re-formed  by  adding  additional  bases  in  that  direction. 

7.1.3  A  Maximum-Likelihood  Spectral  Estimator 

A  non-iterative  solution  for  joint  maximum-likelihood  estimation  of  multiple  independent 
spectral  lines  has  not  been  found.  However,  the  maximum-likelihood  estimate  of  a  single  spectral 
line  can  be  found  conditioned  on  a  previous  estimate  of  the  covariance  matrix  and  any  non-zero 
means.  This  problem  can  be  stated  as  follows:  Given  a  sample  covariance  matrix  S  formed  from 
P  mutually  independent  random  variables 
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V(z*:  -  bkd)(zk  -  6*;d)t  +  y\ 
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(168) 


estimate  a  covariance  matrix  of  the  form  R  -+-  udd^.  R  and  d  are  assumed  to  be  known,  and  a 
is  unknown  but  positive.  The  estimate  of  a  is  found  by  a  maximum-likelihood  estimator.  The 
log-likelihood,  neglecting  constant  terms,  is  proportional  to 


/  =  -  log  R  —  crddt:’  -  tr  (R  -h  ndd*)'^ 


(169) 
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Expanding  the  determinant  as 


;R  4-  rxddf!|  =  j!R!|(l  +  ^R^d) 


(170) 


and  expanding  the  inverse  of  the  covariance  matrix  using  the  Woodbury  identity 


(R  -!-  add1)-1  =  R  1  —  a 


R~1ddtR~1 
1  +  <xdTR_1d 


(171) 


the  terms  of  the  log-likelihood  that  are  functions  of  a  can  be  written 


l  =  -  log(l  —  crdR  Jd)  +  a 


d^'SR-M 

1  -r  adtR  M 


(172) 


Taking  the  derivative  of  this  equation  with  respect  to  a  and  setting  the  result  to  zero  yields 


dt_  -dtRM  _ 

da  1  -f-  crdtR~ xd  l-'-adtR^d 


dfR_1SR_1d  dtR_1SR'1ddtR”1d 

-  a - — — — r~ r^; ~  0 
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and  the  maximizing  spectral  weight  is 


<7=  *  ,y.,  dtR~1SR~1d  -  dtR~1d] 

(dTR  ' 1  d )-  L  J 


(174) 


This  will  also  satisfy  the  necessary  condition  on  the  second  derivative  of  the  log-likelihood.  As  this 
is  assumed  to  be  a  spectral  point,  the  contribution  of  the  energy  must  be  restricted  to  a  positive 
quantity;  therefore,  the  maximum-likelihood  estimate  of  a  which  satisfies  this  constraint  is 


a  =  l  ax(0.  a) 


(175) 


When  the  estimated  spectral  weight  is  non-zero,  then  it  can  be  shown  that  the  increase  in  likelihood 
for  the  covariance  estimate  that  contains  thu  spectral  estimate  compared  to  the  likelihood  without 
this  estimate  is 


A/  =  l( Z  B  R  -  rrddr )  -  l(Z:B.R)  = 


d'R-'SR-'d 
dfR  M 


-  1  -  log 
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(176) 


A /  is  a  monotonic  function  of  (dfR  *SR  ,d)/(d,R  Jd)  since  this  quantity  is  assumed  to  be 
greater  than  or  equal  t^  one  (dtR_1SR'1d  -  d’R'd  >  0)  . 
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This  estimator  can  be  modified  slightly  by  including  an  initial  estimate  of  the  spectral  weight 
in  R.  If  the  initial  estimate  is  v ,  then  Equation  (174)  can  be  written  as  a  function  of  R  and  the 
initial  estimate  so  that 

<t(R)  =d(R  +  zyddt)  +  i/  (177) 

Making  the  restriction  that  a  +  v  be  positive,  the  maximum-likelihood  estimate  of  a  using  the 
initial  estimate  u  is 

a  =  max(0,  a  +  u)  .  (178) 

This  property  will  be  utilized  in  the  proposed  covariance  estimator. 

It  is  interesting  that  the  maximum-likelihood  estimator  for  a  single  point  contains  those  terms 
that  are  used  in  the  update  equation  for  the  spectral  estimates  in  the  EM  iterations  of  Chapter  6. 
The  difference  between  the  update  equation  and  this  estimator  is  that  the  term  added  to  the  current 
estimate  in  the  EM  iteration  is  weighted  by  the  square  of  the  current  spectral  estimate;  while  here, 
the  estimate  is  weighted  bv  the  square  of  the  Capon  spectral  estimator. 

7.2  Description  of  Variable  Basis  Estimator 

In  this  section  an  iterative  algorithm  is  defined  that  can  be  used  to  estimate  the  covariance 
matrices  satisfying  the  spatial  constraint.  The  likelihood  produced  by  the  resulting  sequence  of 
estimates  is  a  non-decreasing  function,  and  stable  points  of  the  estimator  satisfy  the  conditions 
necessary  for  a  maximum-likelihood  estimator. 

This  algorithm  is  based  on  the  estimator  of  Chapter  6  but  does  not  require  the  large  number 
of  terms  that  is  used  to  represent  the  entire  space  of  constrained  covariance  matrices  in  order  to 
represent  the  estimate.  The  data  model  assumed  here  is  the  same  as  that  used  in  the  previous 
chapters.  P  mutually  independent  data  samples  are  available;  the  first  G  are  assumed  to  be  non¬ 
zero  mean  primary  data  z\  . . .  z g.  with  density  Ar(6^d,  R),  and  the  last  K  are  the  secondary  data 
zg^i  •  •  •  zp.  with  density  A:(0,  R). 

The  estimation  algorithm  consists  of  the  five  steps  that  follow. 

1.  Form  an  initial  estimate  of  the  sample  covariance  matrix  as 

Q  p 

s  =  -jj  51  (z*  -  M)(z*  -  M)f  +  ]T  Zkzk  •  (179) 

1  [k=\  k=G+ 1 

where  bk  =  cEz^/d'd.  These  are  the  estimates  of  b  that  result  from  the  assumption 
that  the  noise  covariance  matrix  is  a  diagonal  matrix.  The  term  of  the  sample 
covariance  matrix  corresponding  to  the  primary  data  can  also  be  written  as 


(180) 


Dz*  -  6*dKz*  -  6*d)f  =  (*  -  ^)  E  z*zl  (*  - 

It  will  be  shown  in  the  following  chapter  that  the  terms  in  Equation  (179)  are  suffi¬ 
cient  statistics  for  the  estimate  of  the  covariance  matrix. 


2.  Form  an  initial  set  of  L  basis  functions  and  weights.  It  is  assumed  that  no  prior 
knowledge  of  the  interference  directions  exists.  Simulations  have  shown  that  the 
number  of  iterations  required  to  achieve  a  given  likelihood  is  dependent  upon  the 
initial  estimate:  therefore,  initialization  is  discussed  in  detail  here.  Two  methods  are 
proposed  for  this  initialization.  For  both  of  these  methods,  if  the  constraint  cone 
has  a  hvperplane  boundary  as  shown  in  Chapter  5,  then  the  bases  for  this  boundary 
should  be  in  the  initial  set  of  bases.  Also,  if  the  variance  of  the  receiver  noise  is 
unknown  and  must  be  estimated,  the  elementary  vectors  will  be  used  to  form  N  of 
these  basis  functions. 


•  The  first  method  begins  with  L  bases  that  are  parameterized  by  a  uniform 
distribution  of  angles  throughout  space.  The  direction  vectors  corresponding 
to  these  basis  functions  are  the  columns  of  the  matrix  D.  which  should  be 
full  rank.  The  spectral  weights  on  these  direction  vectors  are  set  to  a  value 
greater  than  zero  but  smaller  than  the  expected  contribution  of  energy  due 
to  receiver  noise  or  due  to  spatially  distributed  noise  for  the  region  near  that 
direction.  With  these  spectral  weights  and  basis  functions,  form  an  estimate  of 
the  covariance  matrix. 


Using  the  set  of  basis  functions,  search  for  the  basis  function  whose  weight,  if 
increased,  would  result  in  the  largest  increase  in  likelihood.  This  is  the  basis 
function  generated  by  the  vector  that  maximizes  Equation  (176)  or  equivalently 


dk  = 


argmax 

dtD 


dtR-1SR-1d 

dlR_1d 


(181) 


For  this  basis  function  and  current  weight  ck.  calculate  the  new  spectral  weight 
a,  =  max  (ak.ak  +  - -  [d[R-1SR-1d*  -  d[R_1djt]  j  .(182) 

V  (d^R  dk)  1  / 

Using  this  spectral  weight,  update  the  inverse  of  the  covariance  matrix  using 
the  Woodbury  formula,  an  0{N2)  operation,  and  update  the  term  R-1SR-1. 
Using  the  Woodbury  formula,  this  can  also  be  implemented  as  an  0(A'2)  op¬ 
eration.  Repeat  this  estimation  of  initial  spectral  weights  until  L  weights  have 
been  set  or  until  the  estimated  weight  for  the  basis  function  with  the  maximum 
increase  in  likelihood  results  in  a  decrease  or  no  change  in  its  initial  weight. 
The  same  weight  will  not  be  updated  consecutively.  The  estimate  is  linear 
with  respect  to  the  spectral  weight,  and  the  gradient  has  been  set  to  zero  by 
the  estimated  weight. 


•  The  second  initialization  method  begins  similarly  to  the  first,  by  initializing 
part  of  the  directions  of  arrival  of  the  spectral  points  uniformly  throughout 
space  with  a  small  weight;  however,  it  then  departs  by  allowing  the  directions 
for  the  other  bases  to  vary.  Using  L/2  fixed  vectors  and  their  initial  weights, 
form  an  estimate  of  the  covariance  matrix  R.  With  this  initial  estimate  of  the 
covariance,  the  remaining  directions  of  arrival  and  weights  will  be  found  by  us¬ 
ing  the  conditional  maximum-likelihood  estimate  for  the  weights.  Rather  than 
maximizing  the  initial  likelihood  using  the  fixed  bases,  the  directions  of  arrival 
for  the  other  1/2  basis  functions  are  varied  so  as  to  maximize  the  likelihood. 


Beginning  with  the  initial  estimate  of  R,  perform  a  search  for  the  direction  of 
arrival  that  would  result  in  the  maximum  increase  in  likelihood.  This  is  the 
direction  (#k,0fc)  that  satisfies  Equation  (176)  or 


(0Ar.d>/c)  = 


argmax 

e,o 


dt(^o)R-1SR-1d(^,o) 

dt(0,o)R-1d(0.0) 


(183) 


This  is  a  computationally  intensive  task  requiring  on  the  order  of  iV2  com¬ 
putations  per  point  in  the  search.  This  can  be  performed  by  using  an  array 
of  fixed-basis  vectors  similar  to  the  array  used  in  the  estimator  of  Chapter  6. 
After  the  direction  of  arrival  for  the  maximum  increase  in  likelihood  is  found 
by  this  coarse  search,  the  direction  can  be  refined  by  the  use  of  a  search  be¬ 
tween  the  directions  of  arrival  of  the  array  of  basis  vectors  until  a  suitable 
stopping  criterion,  such  as  the  difference  in  the  angle  of  arrival  is  within  some 
f.  is  reached.  The  weight  for  this  direction  can  be  found  by  the  conditional 
maximum-likelihood  estimate  as 

&  —  max(0.  — t — - -  [dtR-1SR_1dt  -  dtR-1d*l)  .  (184) 

(d[R-1dfc)2  k  k 

After  estimating  a  direction  of  arrival  and  weight,  the  inverse  of  the  initial 
co%"ariance  estimate  and  the  matrix  product  R_1SR_1  are  updated,  and  the 
process  is  continued  until  all  L  directions  and  weights  have  been  found.  If 
during  this  process  the  weight  for  the  direction  vector  corresponding  to  the 
maximum  increase  in  likelihoc  d  is  negative,  the  remaining  direction  vectors 
should  be  distributed  as  uniformly  throughout  space  as  possible  and  the  weights 
set  to  a  small  value. 


For  both  of  these  initialization  methods,  the  initial  set  of  bases  is  chosen  so  as  to  form 
a  linearly  independent  set.  The  bases  for  all  iterations  are  restricted  to  form  a  linearly 
independent  set.  The  resulting  set  of  hases  forms  an  L-dimensional  constraint  cone. 
With  this  restriction,  there  will  not  be  two  bases  parameterized  by  the  same  direction 
of  arrival.  This  allows  the  estimate  to  move  in  any  direction  interior  to  the  const  raint 
region  during  the  estimation  procedure. 
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3.  Using  this  set  of  basis  functions  and  initial  weights,  perform  EM  iterations  using  the 
method  of  Chapter  6.  These  equations  are  repeated  here: 


£Pjl  =  Sp  +  (Rp  Vrp  1  -  RP  DEP 
Rp+1  -  DSp+1Dt 
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and 


K  = 


dj  R^z, 
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where  the  sample  covariance  matrix  is 
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(188) 


The  direction  vector  corresponding  to  the  mean  does  not  need  to  be  one  of  the  basis 
vectors  that  is  used  to  estimate  the  covariance  matrix  because  the  modification  to 
the  estimator  of  Chapter  6  maximizes  the  likelihood  in  each  iteration  separately  with 
respect  to  the  mean. 


4.  If  the  estimate  is  moving  outside  of  the  polyhedral  cone  formed  by  the  initial  set 
of  basis  functions,  the  estimate  will  approach  the  boundary  of  the  constraint  set. 
This  is  determined  by  examining  the  conditional  maximum-likelihood  spectral  es¬ 
timates  given  the  current  estimate  of  the  covariance  matrix.  If  the  conditional 
maximum-likelihood  estimate  would  be  zero  using  the  modified  form  for  the  con¬ 
ditional  maximum- likelihood  estimate  found  in  Equation  (178),  then  the  estimate  is 
approaching  a  hyperplane  boundary.  Examining  the  approach  of  the  estimate  to  a 
hvperplane  boundary  uses  only  those  quantities  that  need  to  be  calculated  to  update 
the  spectrum  estimates  using  the  EM  algorithm  and  requires  one  additional  multiply 
and  a  divide  per  point.  Examining  the  approach  to  a  hyperplane  boundary  should  be 
accomplished  before  the  update  of  the  spectrum  estimates.  If  any  conditional  spec¬ 
trum  estimate  would  be  zero,  then  the  basis  function  and  weight  should  be  removed 
and  replaced  with  a  basis  function  that  when  added  would  result  in  the  maximum 
increase  in  likelihood  using  the  procedure  outlined  in  the  initialization,  that  is.  with 
the  basis  function  corresponding  to  the  angle  of  arrival  that  maximizes 


(6.  o)  = 


argmax 

9,<t> 


dt(e.<p)R-1SR-1d(0.o) 

d^.tMR-'d^.o) 


(189) 


If  there  is  not  a  different  basis  function  that  would  result  in  an  increase  in  likelihood, 
then  the  current  basis  function  and  weight  are  retained.  It  has  been  found  that 
replacing  a  single  basis  function  and  weight  in  any  iteration  is  sufficient.  Additionally, 
it  has  been  found  that  the  examination  of  the  approach  to  a  hvperplane  boundary 
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does  not  need  to  be  performed  for  each  iteration  since  the  estimate  will  be  interior 
to  the  new  polyhedral  cone  and  can  move  within  this  cone. 

5.  Perform  steps  3  and  4  until  a  suitable  stopping  criterion  is  reached.  The  stopping  cri¬ 
terion  could  be  based  on  the  change  in  likelihood,  or  it  could  be  based  on  performing 
a  fixed  number  of  iterations. 

7.3  Estimator  Properties 

In  Chapter  6  it  was  shown  that  under  certain  conditions  a  maximum-likelihood  estimator 
for  the  mean  and  the  covariance  may  not  exist  for  this  model.  There  will  typically  be  a  maximal 
proper  subspace  orthogonal  to  the  direction  of  arrival  for  the  mean,  as  discussed  in  Section  6.3.  If 
all  of  the  data  are  assumed  to  have  a  non-zero  mean  and  the  sample  covariance  matrix  is  not  full 
rank,  the  likelihood  will  be  unbounded  above  and  the  maximum-likelihood  estimate  will  not  exist. 

For  a  non-singular  sample  Covariance  matrix  and  data  that  are  assumed  to  be  zero  mean,  the 
maximum-likelihood  estimate  will  exist  [28].  If  the  sample  covariance  matrix  is  singular,  then  the 
existence  of  the  maximum-likelihood  estimate  will  be  determined  by  the  array  geometry,  and  the 
dimension  of  the  constraint  space  and  general  statements  will  not  be  made  here. 

To  ensure  that  a  maximum-likelihood  estimate  exists,  receiver  noise  with  a  non-zero  variance 
should  be  included  in  the  model.  This  adds  to  the  covariance  matrix  during  each  iteration  just  as 
it  did  in  Chapter  6. 

A  stable  point  of  this  algorithm  satisfies  the  necessary  conditions  for  a  maximizer  of  the 
likelihood.  If  the  necessary  conditions  on  the  gradient  are  satisfied,  then  the  sequence  of  estimates 
produced  by  the  method  of  Chapter  6  wrill  be  stable  as  the  gradient  of  the  likelihood  appears  directly 
in  the  update  equation.  This  was  pointed  out  in  Chapter  6.  If  these  estimates  are  stable,  then  the 
estimate  is  interior  to  the  finite  bases  constraint  cone,  and  no  additional  bases  will  be  added.  If  the 
gradient  is  positive  in  any  direction  interior  to  the  constraint  set.  then  the  sequence  of  estimates  will 
move  in  that  direction  until  the  estimate  is  stable  or  until  the  estimate  approaches  a  hvperplane 
and  an  additional  basis  dyad  is  added.  Since  a  new  basis  element,  when  added,  maximizes  the 
likelihood,  the  estimate  will  only  move  toward  a  maximizer  of  the  likelihood. 

7.4  Comparisons  to  Fixed-Bases  Estimator 

7.4.1  Convergence 

The  estimation  procedures  presented  in  this  chapter  appear  to  be  much  more  complex  than 
the  fixed-bases  estimation  procedure.  In  this  section  we  show  that,  although  the  variable  basis 
estimator  is  theoretically  more  complex,  the  computational  requirements  can  be  fewer  than  the 
requirements  for  the  fixed-bases  estimator. 

The  typical  environment  for  adaptive  beamforming  and  detection  can  have  a  few  spatially 
discrete  interference  sources  with  energy  levels  that  are  50-80  dB  higher  than  the  receiver  noise 
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[51].  With  levels  this  high,  the  true  covariance  matrix  will  lie  near  the  boundary  of  the  constraint 
region.  In  order  to  adequately  represent  these  matrices,  the  fixed-bases  estimator  must  use  a  large 
number  of  terms  in  the  finite  sum  and  accept  a  small  loss  in  detection  and  beamformer  performance 
or  use  a  small  number  of  terms  and  accept  a  larger  loss  in  performance.  This  performance  trade-off 
is  not  required  if  the  basis  functions  are  allowed  to  vary. 

Reducing  the  number  of  terms  required  to  represent  the  estimate  can  result  in  a  more  rapid 
approach  of  the  estimate  to  the  maximum-likelihood  solution.  Recalling  the  update  equations  for 
the  spectral  estimates  of  Equation  (167) 

of 1  ^(l+o^A)  (190) 

for  a  particular  covariance  matrix,  if  the  energy  is  spread  among  a  larger  number  of  spectral  weights, 
then  the  estimate  will  not  move  as  far  in  the  direction  of  the  gradient.  Assume  that  the  covariance 
estimate  has  an  error  term  Rp  that  is  approaching  zero.  If  this  term  is  formed  by  a  single  spectral 
weight,  then 


RP  =  <7pddf  (191) 

If  this  term  is  represented  by  a  large  number  of  spectral  weights,  then  we  have 

G 

Rf  =  y<dsdJ  (192) 

g~  l 

It  will  be  assumed  that  the  array  geometry  is  such  that  Y^=  i  ag  =  aP-  aru*  tie  spectral  weights  are 
evenly  distributed.  Rp  is  nearly  the  same  in  either  case.  This  could  occur  if  the  direction  vectors 
forming  the  basis  functions  are  closely  spaced.  Then  <tp  =  ap/G.  It  will  also  be  assumed  that  the 
gradient  of  the  likelihood  is  equal  for  all  of  the  terms.  Then 

HT1  =  i:^(l-^A^)dgd],  (193) 

g=l 

for  the  estimate  containing  many  terms,  and 

Rp^'  =  <rp(l  ■+■  A<7p)d9d*  (194) 

for  the  estimate  with  a  single  term.  The  fractional  change  in  the  first  case  is  Acrp/G.  and  for  the 
second  it  is  Acrp,  showing  that  with  the  energy  spread  among  many  directions  will  result  in  slowing 


85 


the  rate  of  convergence.  Re  is  approaching  zero,  so  the  estimate  cannot  go  too  far  in  the  direction 
of  the  gradient  in  either  case. 

7.4.2  Computational  Complexity 

The  algorithm  proposed  in  this  chapter  has  lower  computational  requirements  than  the  fixed- 
bases  estimator.  The  computational  burden  will  be  investigated  in  terms  of  the  number  of  multi¬ 
plications  required.  The  number  of  terms  in  the  fixed-bases  finite  sum  is  a  quantity  that  can  be 
traded  off  against  the  performance  of  detection  and  beamforming  algorithms  using  the  estimates 
produced.  The  number  of  bases  is  a  multiple  of  the  dimension  of  the  constraint  space,  for  example. 
Q.  The  differences  in  the  number  of  multiplications  will  occur  for  the  spectral  weight  updates  and 
generating  the  covariance  from  the  spectral  weights.  The  variable  basis  method  will  also  have  the 
overhead  of  determining  the  weights  when  adding  or  reducing  the  number  of  terms  in  the  spec¬ 
trum.  The  quadratic  term  <!*(•• -)d  of  the  spectrum  update  requires  X2  +  O(N)  multiplications 
for  each  spectral  point  in  each  iteration.  For  the  fixed-bases  esiimator,  this  will  be  QLX 2  +  O(N). 
The  number  of  multiplications  for  the  method  discussed  in  this  chapter  for  the  spectral  update  is 
LN2  +  0(N). 

To  calculate  a  new  covariance  matrix  from  the  spectrum  samples  requires  X2 ,'2  +  O(X) 
multiplies  for  each  spectral  point.  For  the  fixed-bases  estimator  the  number  per  iteration  will  then 
be  QLX2/ 2  +  0( X),  and  for  the  variable  basis.  LX2  2  +  O(X). 

Each  time  a  basis  function  is  added  or  subtracted,  a  search  for  the  maximum  increase  in  like¬ 
lihood  is  performed.  The  most  intensive  task  is  finding  the  approximate  location  of  the  maximum, 
which,  if  performed  using  fixed  bases,  would  require  approximately  QLX 2  —  O(.V)  multiplies.  Af¬ 
ter  this,  increasing  the  resolution  using  a  binary  search  may  require  an  additional  10A'2  multiplies. 
Updating  the  inverse  of  the  covariance  matrix  and  gradient  of  the  likelihood  requires  3X2  opera¬ 
tions  plus  3iV2  for  the  removal  of  the  weight  that  is  approaching  zero.  Additionally,  the  updates 
for  the  spectral  estimates  must  be  recalculated  requiring  LX~  0{ X)  multiplies.  It  is  assumed 
that  a  basis  function  must  be  added  every  J  iterations:  therefore,  the  number  of  multiplications 
per  iteration  is  [QL  —  L  +  16) AT2/ J  +  O(X).  For  each  method,  the  number  of  real  multiplications 
per  iteration,  ignoring  lower  order  terms,  is 


fixed  basis  6QLX2 

variable  basis  6LX2  +  4( QL  +  L  +  16).Y2/J 

Looking  at  the  relative  values  of  these  quantities,  the  terms  Q.  X  and  J  are  all  about  the  same 
magnitude.  Assuming  that  QL  >>  £  +  16,  the  fixed-bases  estimator  will  require  0(LN3)  multiplies, 
while  the  variable  basis  method  will  require  0(LX2).  The  variable  basis  method  actually  requires 
less  computation  than  the  fixed-bases  method. 
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7.5  Simulation  Results 


A  few  brief  simulations  will  be  used  to  compare  the  estimator  of  this  chapter  to  the  fixed-bases 
estimator  of  Chapter  5.  A  4-element  minimum  redundancy  linear  array  has  been  simulated.  The 
dimension  of  the  space  of  matrices  is  13.  The  fixed-bases  estimator  uses  26,  130,  and  900  direction 
vectors  to  form  the  basis.  All  of  the  estimators  also  estimated  a  diagonal  term  that  was  constrained 
to  be  identically  distributed  for  each  of  the  vector  elements.  The  variable  basis  estimator  of  this 
chapter  used  a  fixed  resolution  of  0.05  degree  or  the  equivalent  of  3600  basis  vectors.  The  true 
covariance  was  used  in  place  of  the  sample  covariance  so  that  the  maximum  log-likelihood  was 
known,  and  this  was  subtracted  from  the  likelihoods  before  plotting  so  that  the  maximum  for  each 
of  the  plots  is  0.  The  need  to  adjust  the  bases  was  determined  every  fourth  iteration  for  the  method 
of  this  chapter.  Two  interference  environments  are  simulated.  The  first  environment  consists  of 
receiver  noise  with  a  variance  of  1,  uniform  spatial  interference  with  a  variance  of  1  for  a  single 
vector  element  (which  contributes  an  identity  term  here  due  to  simulated  sensor  placement),  and 
three  discrete  interference  sources  located  at  azimuths  of  -25.78.  16.04.  and  34.95  degrees  with 
power  levels  of  30.  40.  and  60  dB  across  the  array,  respectively,  when  compared  to  the  receiver 
noise.  The  second  environment  consists  only  of  the  receiver  noise  and  uniform  spatial  interference, 
and  the  covariance  matrix  is  interior  to  the  polyhedral  cone  formed  by  the  set  of  bases  for  each  of 
the  fixed-bases  estimators. 

A  plot  of  the  likelihood  versus  iteration  is  shown  here  in  Figure  17.  We  see  that  the  variable 
bases  estimator  has  a  more  rapid  rate  of  convergence  and  results  in  a  final  likelihood  that  is  greater 
than  the  fixed-bases  estimator.  Iterations  where  basis  functions  were  added  are  shown  by  the 
increase  in  the  rate  of  convergence.  The  estimator  of  this  chapter  re-formed  the  set  of  bases  15 
times  during  these  iterations  for  the  first  interference  set  and  did  not  require  changing  the  bases 
for  the  second.  The  estimators  of  Chapter  6  were  initialized  by  a  constant  spectrum  value  (all 
ones).  The  difference  in  the  convergence  rate  is  quite  dramatic  when  comparing  the  method  of  this 
chapter  to  the  estimator  of  the  previous  chapter. 

To  show  that  the  performance  difference  is  not  due  solely  to  the  new  initialization  method,  the 
estimator  of  Chapter  6  was  initialized  by  the  first  method  described  in  this  chapter  for  Figure  18. 
Although  this  initialization  results  in  a  higher  initial  likelihood,  the  estimator  with  the  large  number 
of  terms  (with  some  spectral  weights  approaching  zero  or  with  the  small  number  of  terms  that 
cannot  represent  the  true  covariance  matrix)  does  not  converge  after  1000  iterations.  One  additional 
plot.  Figure  19.  will  be  used  simply  to  expand  the  scale  of  the  second  plot  in  Figure  18.  This  shows 
that  the  initialization  causes  the  estimate  to  be  very  nearly  the  maximum-likelihood  estimate  for 
this  noise  environment  and  the  small  number  of  bases  used. 

For  the  second  interference  environment,  all  of  the  fixed-bases  estimators  are  able  to  represent 
the  covariance  matrix:  however,  for  the  estimator  with  900  bases  the  convergence  is  so  slow  that 
the  estimator  has  not  converged  after  1000  iterations. 

These  plots  show  the  trade-off  between  the  number  of  bases  of  the  previous  chapter  and  the 
ability  to  represent  strong  discrete  interference  sources  with  the  rate  of  convergence. 
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Likelihood  vs  iteration 


iteration 


Figure  19.  Comparison  of  convergence  rates  expanded  scale. 


7.6  Conclusion 

In  Chapter  5  the  maximal  representation  of  covariance  matrices  was  discussed.  In  this  chapter, 
an  estimation  method  has  been  presented  that  is  based  on  the  maximal  representation.  This  method 
has  been  shown  to  have  computational  and  convergence  benefits  compared  to  the  fixed-basis  finite 
representation  estimator  developed  in  Chapter  6. 

In  the  following  chapter,  we  will  use  the  covariance  estimates  and  the  likelihood  found  through 
this  procedure  to  develop  adaptive  detection  and  beamforming  methods. 
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8.  BtAMFORMING  AND  DETECTION  METHODS 


The  material  presented  in  the  three  previous  chapters  focused  on  estimating  the  parameters 
that  are  indicated  by  the  signal  model  without  explicitly  stating  how  the  results  of  these  estimation 
algorithms  would  be  used.  In  this  chapter,  several  beamforming  and  detection  methods  that  make 
use  of  these  parameter  estimates  and  estimation  algorithms  are  presented. 


8.1  Beamforming 


In  this  section,  methods  of  estimating  the  mean  corresponding  to  a  given  direction  are  pro¬ 
posed.  These  beamforming  methods  are  based  on  the  approaches  discussed  in  Chapter  3  and  use 
maximum-likelihood  techniques.  It  is  assumed  that  a  single  primary  vector  z  may  have  a  non-zero 
mean  of  the  form  6d.  and  K  mutually  independent  secondary  vectors  are  available  that  are  subject 
to  a  density  function  with  the  same  covariance  as  the  primary  vector.  The  object  of  the  beamformer 
is  to  estimate  the  complex  amplitude  b. 


1.  The  first  method  is  to  perform  joint  estimation  of  the  mean  and  the  covariance  matrix 
by  the  maximum-likelihood  estimation  procedures  introduced  in  Chapters  6  and  7. 
These  procedures  will  use  all  of  the  available  data  to  form  the  estimates  of  the  mean 
and  the  covariance. 


2.  The  second  method  is  similar  to  that  suggested  by  Reed  et  al.  { 10! :  Substitute 
the  constrained  estimate  of  the  covariance  using  only  the  secondary  data  for  the 
covariance  matrix  in  the  maximum  signai-to- noise  ratio  beamformer 

w  =  kRld  (195) 


Additionally,  the  constant  k  is  set  so  that  the  gain  in  the  desired  look  direction  is 
equal  to  unity.  The  weight  vector  is  then 

(.96) 


w 


d '  R  ‘d 


This  weight  vector  will  be  used  to  form  an  estimate  of  the  mean  in  the  primary  vector 
z  by 


b  =  w*z 


(197) 


For  unconstrained  maximum-likelihood  covariance  estimates,  these  two  procedures  produce  equiva¬ 
lent  estimates  of  the  mean.  This  property  has  not  been  proven  to  be  true  for  constrained  covariance 
estimates. 


Beamformers  based  on  constrained  covariance  estimates  will  be  compared  to  the  beamformer 
using  unconstrained  maximum-likelihood  covariance  estimates  and  to  the  optimum  linear  beam- 
former  formed  when  the  covariance  matrix  is  known.  Because  closed-form  expressions  for  the 
covariance  estimates  and  hence  the  mean  estimates  have  not  been  found,  the  beamforming  methods 
will  be  compared  through  simulations. 


Performance  will  be  evaluated  by  comparing  the  biar  and  variance  of  the  estimates  produced 
by  each  method.  This  also  compares  the  variance  of  the  estimates  to  the  Cramer-Rao  bound  when 
the  co%"ariance  is  known  since  the  known  covariance  maximum-likelihood  estimator  is  efficient. 

The  covariance  estimates  will  be  evaluated  in  terms  of  their  application  to  heamforming  by 
comparing  the  signal-to- noise  ratio  loss  factor  introduced  by  Reed  et  al.  [10]  and  discussed  in 
Chapter  3.  The  loss  factor  is 

IwM'2 

P  =  cPR-MwTRw 


I^R-^!2 

cPR-MdtR-iRR-M 


(198) 


This  loss  factor  will  verify  the  comparison  of  the  variance  since,  for  an  unbiased  estimator,  minimiz¬ 
ing  the  variance  is  equivalent  to  maximizing  the  sigr al-to-noise  ratio.  The  loss  factor  is  a  random 
variable  and  will  be  presented  as  a  histogram  to  indicate  the  corresponding  density  function. 

Additionally,  the  spatial  response  of  the  beamforming  methods  will  be  compared.  This  is 
presented  simply  to  show’  the  variance  in  the  spatial  response  produced  by  each  of  the  beamforming 
methods. 


8.2  Adaptive  Detection 


The  proposed  detectors  assume  that  there  is  a  single  primary  vector  z.  which  may  have  a 
non-zero  mean  (G=l)  of  the  form  6d.  d  is  known,  and  b  is  assumed  to  be  unknown.  K  additional 
mutually  independent  secondary  vectors  Z  =  {zi  ■  ■  -  Zk}  are  also  available,  and  these  vectors  share 
the  same  covariance  matrix  R. 


Two  methods  of  utilizing  the  algorithms  presented  in  Chapters  6  and  7  are  introduced. 

1.  The  first  detection  method  is  to  utilize  the  generalized  likelihood  ratio  test  procedure 
[2]  given  by 


A(z,  Z) 


max  /( z,  Z:  b ,  R\Hi) 

6.R!//, 


max  /( z,  Z;  Ri  H0) 
R\H0 


lh 

Ho 


(199) 


The  test  statistic  for  this  method  cannot  be  written  in  closed  form;  rather,  the 
likelihood  ratio  must  be  computed  explicitly  from  the  maximum-likelihood  estimation 
procedures  under  each  hypothesis. 


2.  The  constrained  maximum-likelihood  estimate  of  the  covariance  matrix  based  on 
the  secondary  data  may  be  substituted  for  the  known  covariance  in  the  unknown 
mean,  known  covariance  generalized  likelihood  ratio  test.  This  can  be  justified  if  the 
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estimate  of  the  covariance  matrix  is  estimated  with  sufficient  accuracy.  The  resulting 
test  is 


dTR  1  z 
dtR_1d 


2  Hi 

~  ^  1 


Ho 


(200) 


For  both  of  these  methods,  the  test  statistics  are  positive.  In  the  simulations  discussed  in  the 
following  chapter,  a  monotonic  function  of  the  test  statistics  is  plotted. 

Without  a  closed-form  expression  for  the  test  statistics  under  these  methods  rWi-w  rUncjty 
functions  and  evaluating  performance  based  on  these  density  functions  cannot  be  accomplished. 
However,  some  of  the  properties  of  these  detectors  will  be  discussed  in  the  next  section.  How  these 
detectors  perform  compared  to  the  other  detectors  will  be  evaluated  by  simulations  in  the  next 
chapter. 


8.3  Discussion 

8.3.1  Beamforming 

Because  closed-form  expressions  for  the  constrained  covariance  estimates  have  not  been  found, 
an  explicit  formula  for  the  loss  in  signal-to-noise  ratio  has  not  been  derived.  Other  properties  of  the 
beamformers  can  be  proven.  For  the  first  beamformer.  the  sufficient  statistic  for  the  estimate  of  the 
covariance  matrix  only  involves  the  secondary  data  and  the  primary  data  that  are  orthogonal  to 
the  subspace  spanned  by  the  known  direction  vector.  In  addition,  in  the  vector  space  of  matrices, 
the  primary  data  will  only  affect  the  estimated  covariance  matrix  in  the  subspace  orthogonal  to 
ddF  It  can  be  shown  that  the  second  beamformer  provides  an  unbiased  estimate  of  the  mean. 

The  independence  of  the  estimated  covariance  matrix  to  that  portion  of  the  primary  data  in 
the  span  of  the  direction  vector  can  be  shown  by  writing  the  log-likelihood  and  then  maximizing 
with  respect  to  the  amplitude  of  the  mean  b.  Substituting  the  maximum-likelihood  estimate  of  b , 
the  log-likelihood  without  simplifying  is 


N P  log  7T  -  P  log  ijRjj 
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lZfc  .(201) 


The  terms  corresponding  to  data  vectors  that  are  non-zero  mean  can  be  written 
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dfR-1z*  i  (  dtR-^ 

Zfc  -  dtR^Id  dJ  R  ^  “  dtR^d  d 


=  zlR  i  I- 


R-iddfR"i 

dtR-M 


R-*ddfR-*\  i 

•-“dlTTd-  R‘,z* 


=  ^R"i(i-Rdmd-id"i)R^Zt  -  (202) 

The  terms  in  the  braces  (•)  on  the  right  are  projection  matrices  with  a  null  space  for  components 
in  the  primary  vector  in  the  subspace  spanned  by  d.  This  can  be  clearly  seen  by  the  equality 


slR-i  I 


R-^ddfR-i 

dtRM 


R*  I 


R“  4(1- 


R  iddfR  ^ 
dtR-M 


The  result  is  that  the  maximum-likelihood  estimate  of  the  covariance  matrix  only  depends  on  the 
vector  components  of  the  non-zero  mean  data  in  the  subspace  orthogonal  to  d.  From  this  it  can  be 
seen  that  the  sufficient  statistic  for  the  estimate  of  the  covariance  matrix  is  the  sum  of  the  outer 
products  of  the  zero  mean  data  and  the  sum  of  the  outer  products  of  the  vector  components  of  the 
non-zero  mean  data  in  the  subspace  orthogonal  to  d. 

A  property  of  the  estimates  is  that,  over  the  vector  space  of  matrices,  the  primary  data  will  not 
affect  the  covariance  estimates  for  variations  of  the  form  <5R  =  crdd^.  The  likelihood  for  variations 
of  this  type  is  not  affected  by  the  primary  data;  the  primary  data  will  only  affect  the  estimate  of 
the  covariance  matrix  for  directions  that  are  orthogonal  to  dd*.  The  necessary  condition  that  must 
be  satisfied  by  maximum-likelihood  covariance  estimates  interior  to  the  constraint  region  is 

61  =  tr  [(R-'SR’1  -  PR-^^r]  - 

S(z‘-i<l^d),R'liRR',(z‘-^l:^d)  =  0  (204) 

k  =  l 

where  S  =  ZZ*.  Now’  let  <5R  =  £crddt.  The  resulting  equation, 


df(R_1SR" 
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does  not  contain  the  primary  data.  Because  the  primary  data  does  not  affect  the  likelihood  for 
variations  of  this  form,  then  this  data  will  not  affect  the  estimate  of  the  covariance  matrix  in  the 
subspace  spanned  by  ddt.  This  can  also  be  shown  by  letting  R  =  R0  +  crddt,  assuming  that  R0 
is  non-singular.  Substituting  this  expression  in  Equation  (201),  using  the  Woodbury  formula  and 
simplifying,  the  log-likelihood  terms  corresponding  to  the  non-zero  mean  data  are  of  the  form 


dtR~1zfc]2 
dtRo"1 


(206) 


a  drops  out  of  the  log-likelihood  expression  for  the  non-zero  mean  terms,  and  the  estimate  of  a 
is  indeterminate  with  respect  to  the  primary  data.  This  was  discussed  in  a  different  context  in 
Chapter  6  where  it  was  shown  that  a  maximum-likelihood  estimator  of  the  covariance  matrix  may 
not  exist  when  all  of  the  data  are  assumed  to  contain  a  mean. 

The  error  in  the  beamformer  plots  used  to  illustrate  the  variance  of  the  estimators  can  be 
related  to  the  covariance  matrix  estimation  error  E.  If  the  true  and  estimated  covariance  matrices 
are  R  and  R.  then 


R  =  R  +  E  (207) 

Since  E  is  Hermitian.  then  it  is  unitarily  similar  to  a  diagonal  matrix  and  has  a  decomposition 

E  =  UEUf  ,  (208) 

where  E  is  diagonal  but  not  necessarily  positive.  The  beamformer  that  results  from  using  this 
estimated  covariance  matrix,  neglecting  the  constant  scale  factor,  is 


w  =  Rld  =  (R  +  UEUf)_1d 


(209) 


Using  the  Woodbury  inversion  formula  [36] ,  the  weight  vector  is 

w  =  R-1d  -  R_1U(UtU  +  E-1)  “1UtR“1d 

=  w0  +  we  ,  (210) 

and  the  magnitude  squared  of  the  array  response  is 

:wfd(i9.  d)'2  =  |(wc  ->r  we)+d {9,  o)|2  =  |w*d(0,  c>)|2  +  w^d(0.  o)d(0,  Ojt(w0  +  we)  .  (211) 
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This  is  the  optimal  response  perturbed  by  the  error  in  the  estimation  of  the  covariance  matrix  that 
appears  in  the  weight  error  vector  we.  If  this  weight  error  vector  is  within  the  span  of  w0,  then  the 
error  will  not  change  the  array  sidelobe  pattern  but  only  the  scaling,  which  would  be  incorporated 
into  the  normalization  of  the  weight  vector. 

The  proof  that  the  second  estimator  is  unbiased  is  based  on  the  statistical  independence  of  the 
primary  and  the  secondary  data  vectors.  The  density  function  of  the  estimated  covariance  matrix 
is  unknown;  however,  the  estimation  of  the  covariance  matrix  does  not  involve  the  primary  data 
vector.  Let  the  mapping  of  the  secondary  data  to  the  estimated  covariance  be  written  as  ©(Z), 
then  the  expected  value  of  the  estimate  is 


df0(Z)z  1  =  [dt0(Z)d] 
dt0(Z)d  J  \dt0(Z)dJ 


and  the  estimator  is  unbiased. 


(212) 


8.3.2  Detection 

Although  closed-form  expressions  cannot  be  found  for  the  estimates  of  the  covariance  matrices, 
some  properties  of  the  detectors  based  on  the  constrained  covariance  estimates  can  be  analyzed.  It 
wiil  be  shown  that  both  of  the  detectors  are  invariant  to  scaling  or  a  change  of  coordinates  of  the 
data  and  are  invariant  to  scaling  of  the  direction  vectors.  Invariance  of  the  test  statistic  to  scaling 
of  the  data  is  necessary  for  a  detector  to  be  CFAR  and  is  itself  a  weak  CFAR  property. 

In  either  of  the  tests  described  earlier,  the  test  statistic  will  be  invariant  to  scaling  of  the  signal 
direction  vector.  This  is  a  desirable  property  since  no  explicit  normalization  has  been  assumed. 
Without  this  property,  detection  performance  would  depend  on  the  scaling.  This  property  can  be 
seen  by  examining  the  test  statistic  for  the  constrained  AMF  test 


!dtR.~1z'2 

cF  R  1  d 


Hi 

£  Q 


Ho 


(213) 


Any  scaling  of  the  direction  vectors  will  be  canceled  in  the  resulting  test  statistic. 

For  the  generalized  likelihood  ratio  test,  this  property  is  based  on  the  form  of  the  log-likelihood 
under  H\.  After  maximizing  the  log-likelihood  with  respect  to  the  mean,  the  terms  in  the  likelihood 
that  arc  dependent  upon  the  signal  direction  vector  are  of  the  form 


cFR-'z 

dtR->dd 


(214) 
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It  can  be  seen  that  the  direction  vectors  enter  into  this  equation  so  that  these  terms  are  invariant 
to  scaling,  and  the  result  is  that  the  maximum  of  the  likelihood  will  be  independent  of  the  scaling 
of  the  direction  vector. 

It  wili  now  be  shown  that  the  test  statistics  are  invariant  with  respect  to  an  invertible  change 
of  coordinates.  An  invertible  change  in  coordinates  will  preserve  the  dimension  of  the  constraint 
space  but  may  not  preserve  the  obvious  structure  of  the  covariance  matrix,  e.g.,  Toeplitz.  Beginning 
with  the  log-likelihood  of  the  numerator, 
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—  NP  log  n  -  P  log  jjRli 
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d)tR-1(zfc  - 


d'R'zk 

dtRM 


d)-  Z  ztR  lz* 

k=G~l 


(215) 


a  necessary  condition  for  the  covariance  matrix,  interior  to  the  constraint  region,  that  maximizes 
the  likelihood  is 
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tr  [(R^SR”1  -  PR~l)6R 
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i>* 
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dTR-1Zfc 
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d^R-1  6RR~\zk 


d'R~lzk 

dfR_1d 


d)  =  0 


(216) 


where  S  =  ^iT=G-i  zk*k-  Based  on  this  expression,  invariance  to  the  scale  and  to  a  change  nf 
coordinates  can  be  shown.  Assume  for  data  set  Z  that  the  solution  to  the  above  equation  is  R. 
Now  subject  the  data  to  an  invertible  transformation  of  coordinates  so  that  yk  =  Vz^,  p  =  Vd. 
For  this  new  data  set  the  covariance  matrix  will  be  a  member  of  the  constraint  space 


Qy  =  {Ry  :  Ry  =  VRzVt,  VR*  G  nx}  (217) 

It  can  be  shown  that  a  solution  to  Equation  (216)  will  be  VRV*  for  data  set  Y.  Substituting  this 
into  Equation  (216),  all  of  the  V  cancel  leaving  the  equation  satisfied  by  R.  This  will  be  true  for 
the  denominator  of  the  generalized  likelihood  ratio  test  as  well.  If  R//0  is  a  maximum-likelihood 
estimate  of  the  covariance  matrix  under  H o.  after  a  change  of  coordinates,  then  VRi^V*  will 
be  a  maximum-likelihood  estimate.  Inserting  these  estimates  into  the  generalized  likelihood  ratio 
formula  leaves  the  arguments  of  the  exponentials  unchanged.  The  determinants  may  be  factored 

IVRV1!  =  |Vj  |R;  [V+|  (218) 
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and  canceled  between  the  numerator  and  denominator,  yielding  the  result  that  the  test  statistic  is 
invariant  with  respect  to  the  transformation  of  coordinates.  One  such  transformation  of  coordinates 
is  a  change  of  scale;  therefore,  the  test  statistic  is  invariant  with  respect  to  a  change  in  scale. 

Using  the  above  proof,  it  can  also  be  shown  that  the  transformation  of  coordinates  for  the 
constrained  adaptive  matched  filter  will  be  canceled  in  the  test  statistic,  and  the  test  statistic  is 
invariant  with  respect  to  the  transformation  of  coordinates. 

Invariance  to  scale  is  a  known  property  of  the  generalized  likelihood  ratio  test  [52].  This 
invariance  to  a  change  in  coordinates  does  not  show  that  the  probability  of  false  alarm  will  be 
independent  of  the  unknown  covariance.  With  invariance  to  the  scale  of  the  direction  vector,  it 
shows  that  the  tests  are  invariant  to  the  scale  but  not  necessarily  the  structure  of  the  true  covariance 
matrix. 


Independence  of  the  density  function  of  the  test  statistic  to  the  structure  of  the  true  covariance 
matrix  is  an  open  research  question.  There  is  an  interesting  property  of  the  test  statistics  that, 
if  it  were  possible  to  condition  on  the  estimated  covariance  matrix,  would  indicate  independence 
of  the  test  statistics  to  the  structure  of  the  true  covariance  matrix.  This  property  holds  when  the 
following  conditions  are  met. 


•  The  sample  covariance  is  not  a  maximum-likelihood  estimate  of  the  covariance.  Since 
the  dimension  of  the  constraint  space  is  less  than  the  dimension  of  the  space  of 
possible  sample  covariance  matrices,  this  typically  holds  with  probability  1. 


•  The  estimate  of  the  constrained  covariance  matrix  is  within  the  interim  of  the  con¬ 
straint  region.  When  this  is  true,  then  the  necessary  condition  for  the  maximizer  of 


•  The  possible  variations  of  the  covariance  matrix  must  include  the  outer  product  of 
the  signal  direction  vector. 


When  the  above  conditions  are  met.  then  Equation  (216)  is  satisfied,  and  variations  of  the  form 
6crddt  can  be  written 


^dtR_1SR_1d  =  d^'d 


(219) 


The  term  on  the  right  appears  in  both  of  the  tests  in  the  argument  of  the  exponent  of  the  likelihood 
equation  under  H\  and  in  the  denominator  of  the  constrained  adaptive  matched  filter.  Substituting 
for  the  denominator  of  the  constrained  adaptive  matched-filter  test  yields 


dtR1zz1R  ‘d  7 
dtR  'SR-M  J  P 

Ho 


(220) 
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If  this  test  could  be  conditioned  on  the  estimate  of  the  covariance  matrix,  it  would  be  easy  to  show 
that  the  resulting  test  statistic  is  independent  of  the  true  and  estimated  covariance  matrices.  As  it 
is  here,  it  serves  to  show  that  the  same  transformation  on  the  primary  and  secondary  data  will  be 
used  to  obtain  the  test  statistic.  This  results  in  a  comparison  of  the  power  of  the  primary  data  in 
the  direction  R-1d  to  the  average  power  of  the  secondary  data  in  the  same  direction. 

In  the  next  chapter  the  methods  presented  in  this  chapter  will  be  compared  to  the  beamform¬ 
ing  and  detection  methods  that  do  not  use  structured  covariance  approaches. 
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9.  ADAPTIVE  BEAMFORMING  AND  DETECTION  RESULTS 


9.1  Introduction 

The  beamforming  and  detection  methods  that  have  been  proposed  in  this  report  and  some  of 
the  methods  discussed  in  Chapter  3  will  be  compared  in  this  chapter.  It  will  be  shown  that  using  the 
knowledge  that  the  true  covariance  matrix  is  structured,  in  the  derivation  of  the  signal-processing 
algorithms,  can  lead  to  a  dramatic  increase  in  performance. 

An  analytic  expression  has  been  found  for  only  one  of  the  structured  covariance  methods 
(Appendix  Ai.  In  order  to  show  that  structured  covariance  methods  may  be  applied  to  adaptive 
beamforming  and  detection,  the  methods  for  which  analytic  expressions  have  not  been  found  will 
be  simulated.  These  simulations  will  be  used  to  compare  the  performance  of  the  adaptive  signal¬ 
processing  algorithms  in  order  to  provide  insight  into  how  these  methods  would  behavt . 

The  parameters  that  may  he  varied  for  the  simulations  are  common  to  both  the  beamformers 
and  the  detector  methods  The  parameters  that  may  be  varied  are 

•  the  number  of  sensors  .V. 

•  the  sensor  locations. 

•  the  interference  environment. 

•  the  number  of  zero  mean  data  vectors  A 

•  the  receiver  noise  level  diagorf  .  .  .  <r\  ). 

•  the  number  of  terms  in  the  finite  sum  approximation  M. 

•  whether  the  constrained  approaches  are  given  knowledge  of  the  receiver  noise  level 
and  structure,  and 

•  the  direction  corresponding  to  the  mean  term. 

Rather  than  attempt  to  illustrate  the  performance  for  all  combinations  of  the  above  param¬ 
eters.  the  simulations  will  be  restricted  to  two  arrays  and  a  fixed  number  of  terms  in  the  finite 
sum.  The  number  of  zero  mean  data  vectors  will  also  be  fixed  at  A  =  A  The  primary  interest 
in  these  simulations  is  to  show  how  the  beamforming  and  detection  methods  perform  in  different 
noise  environments. 

The  first  array  that  will  be  simulated  is  a  minimum-redundancy  linear  array  with  element 
s  fracings 


Array  1  -  (I.  I)  5.  2.0.  3.0 


(221) 
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This  vector  indicates  the  location  of  the  array  elements  in  wavelengths  along  the  array  axis.  The 
covariance  matrix  for  this  array  has  elements  corresponding  to  six  different  non-zero  "lags."  The 
dimension  of  the  constraint  space  is  13.  The  dimension  of  the  entire  space  of  4  by  4  complex 
covariance  matrices  is  16.  The  only  repeated  covariance  elements  are  along  the  diagonal  of  the 
constrained  covariance  matrices.  The  number  of  terms  in  the  finite  sum  approximation  was  fixed 
at  twice  the  dimension  of  the  constraint  space,  and  a  diagonal  term  was  estimated  independently 
for  both  constrained  methods.  The  resolution  of  the  variable  basis  estimator  was  set  to  0.1  degree. 

The  second  array  that  is  simulated  is  a  thinned  linear  array  or  a  uniform  linear  array  with 
elements  removed.  The  element  spacings  are 

Array'2  =  0.  0.5.  1  5.  2.0.  3.0.  3.5.  4.5.  5.0  (222) 


The  dimension  of  the  constraint  space  in  this  case  is  21.  The  resolution  of  the  variable  basis 
estimator  was  set  to  0.05  degree,  and  the  number  of  terms  in  the  fixed  set  of  bases  is  42. 

These  arrays  are  linear  arrays  and  cannot  differentiate  among  sources  that  are  along  vectors 
rotated  around  the  axis  of  the  array;  there  is  a  "cone  of  ambiguity."  Searching  for  directions 
and  bases  that  maximize  the  likelihood  then  reduces  from  the  two-dimensional  search  to  a  one¬ 
dimensional  search. 

Several  different  interference  environments  are  simulated.  These  environments  consist  of 
discrete  sources  plus  a  diffuse  background  interference.  The  diffuse  background  interference  will 
add  a  diagonal  term  to  the  covariance  matrix,  as  will  the  receiver  noise.  The  power  level  and 
location  of  the  discrete  sources  that  have  been  simulated  for  the  4-element  array  are  shown  in 
Table  3.  The  interference  locations  for  the  first  three  environments  are  represented  by  the  basis 


TABLE  3 

Interference  Directions  of  Arrival  and  Intensity;  4-Element  Array 


Model 

Int. 

Sources 

0 

Power 

f) 

Power 

A 

0 

B 

1 

61  8 

.0 

C 

2 

-34  80 

40 

61  80 

50 

D 

1 

57.3 

60 
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vectors  that  are  used  in  the  method  of  Chapter  6  but  not  by  the  initial  bases  for  the  method  of 
Chapter  7.  Environment  D  is  not  within  the  polyhedral  cone  formed  by  the  fixed-bases  estimator 
and  is  not  in  the  initial  cone  formed  by  the  variable  bases  estimator.  For  the  8-element  array, 
slightly  different  environments  were  simulated,  and  these  are  shown  in  Table  4.  The  different 


TABLE  4 


Interference  Directions  of  Arrival  and  Intensity;  8- Element  Array 


Model 

Int. 

Sources 

e 

Power 

e 

Power 

d 

Power 

A 

0 

B 

2 

58.75 

40 

-47.5 

50 

C 

6 

-68.75 

■Hi 

-34.75 

30 

20.5 

30 

37  5 

mm 

41.75 

40 

46.0 

40 

interference  environments  were  used  to  provide  very  different  noise  covariance  matrices  so  that  the 
false  alarm  rate  properties  of  the  detectors  can  be  investigated.  The  number  of  iterations  of  the 
iterative  algorithms  was  fixed  at  100. 

9.2  Beamforming 

The  goal  of  the  beamforming  algorithms  is  to  provide  an  estimate  of  the  mean  b  while  sat¬ 
isfying  some  optimality  criterion.  The  criterion  that  will  be  used  for  comparison  of  the  different 
beamforming  methods  is  the  output  signal-to-noise  ratio.  Four  different  methods  will  be  compared 
in  this  section.  A  summary  of  these  methods  is  given  below. 

•  The  known  noise  covariance  maximum-likelihood  estimate  of  the  mean.  This  method 
is  an  efficient  estimator  of  the  mean.  The  curves  for  this  method  will  be  labeled  ML. 

•  The  unknown  noise  covariance  maximum-likelihood  estimate  of  the  mean.  This 
method  is  equivalent  to  estimating  the  unstructured  covariance  matrix  using  the 
zero  mean  data  samples  and  then  using  this  estimate  in  the  known  noise  covariance 
maximum-likelihood  estimator.  The  curves  for  this  method  will  be  labeled  AML. 

•  Estimate  a  structured  covariance  matrix  satisfying  the  constraints  imposed  by  the 
array  geometry  and  operating  environment,  using  only  the  data  that  is  known  to  be 
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zero  mean.  Use  this  estimate  in  the  known  covariance  maximum-likelihood  estimator 
of  the  mean.  The  curves  here  will  be  labeled  CML6  and  CML7  for  the  estimation 
techniques  of  Chapters  6  and  7,  respectively. 

•  Estimate  jointly  the  structured  covariance  matrix  and  the  mean  using  all  of  the  data 
in  a  single  estimation  algorithm.  The  curves  here  will  be  labeled  J6  and  J7  for  the 
methods  of  Chapters  6  and  7,  respectively. 

The  covariance  matrix  estimation  procedure  based  on  the  one-dimensional  constraint  space  will  not 
be  discussed  as  a  beamforming  method.  The  signal-to-noise  ratio  will  not  exhibit  a  loss  for  this 
method  since,  assuming  that  the  structure  is  correct,  the  signal-to-noise  ratio  is  independent  of  the 
estimated  covariance  matrix.  This  is  discussed  in  Appendix  A. 

9.2.1  Spatial  Response 

The  first  three  beamforming  methods  provide  estimates  of  the  mean  by  a  linear  filtering  op¬ 
eration  on  the  non-zero  mean  data  vector.  Comparing  the  spatial  response  of  these  beamformers 
provides  insight  into  how  well  the  covariance  estimate  has  been  estimated.  The  known  covariance 
maximum  signal-to-noise-ratio  linear  filter  response  is  shown  in  Figure  20  for  environment  C.  The 
responses  for  six  realizations  of  the  linear  filter  using  the  unconstrained  maximum-likelihood  co- 
variance  estimate  are  shown  in  Figure  21.  Figures  22  and  23  show  six  realizations  of  the  filter 
using  constrained  maximum-likelihood  covariance  estimate  of  Chapters  6  and  7.  respectively.  It 
should  be  emphasized  that  the  sidelobes  do  not  reflect  the  loss  in  signal-to-noise  ratio  that  might 
be  associated  with  a  non-adaptive  system  for  the  same  level  of  sidelobes.  We  can  get  a  sense  of  the 
loss  in  signal-to-noise  ratio  by  comparing  the  adaptive  beamformers  to  the  optimal:  however,  high 
sidelobes  in  themselves  do  not  indicate  poor  performance. 

Figures  24,  25,  26,  and  27  show  the  beamformer  response  for  the  8-element  array  and  inter¬ 
ference  environment  C  where  the  difference  in  the  beamformer  responses  is  more  obvious.  In  the 
next  section  the  loss  factor  corresponding  to  these  beamformer  responses  will  be  shown. 

9.2.2  Signal-to-Noise  Ratio  Loss  Factor 

In  this  section  the  loss  in  signal-to-noise  ratio  that  would  occur  through  the  use  of  the  covari¬ 
ance  estimates  rather  than  the  true  covariance  matrix  in  a  linear  beamformer  will  be  compared. 
The  density  of  this  loss  is  known  for  the  unconstrained  covariance  estimates,  so  the  histogram  of 
200  realizations  of  the  constrained  covariance  methods  will  be  compared  to  this  density  function. 
Figures  28  and  29  are  histograms  of  the  loss  factor  from  the  simulations  of  the  last  two  methods 
overlaid  with  the  beta  density  function,  which  is  the  loss  factor  density  when  the  beamformer  is 
based  on  the  unstructured  covariance  estimate.  The  histograms  show  that  the  constrained  covari¬ 
ance  estimates  exhibit  lower  loss  than  the  beamformer  based  on  the  unconstrained  estimates.  In 
these  histograms,  the  true  covariance  was  in  the  polyhedral  cone  formed  by  the  finite  basis  used  to 
estimate  the  structured  covariance. 
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Figure  26.  Beamformer  response  us¬ 
ing  Chapter  6  structured  covariance  es- 


Figure  27.  Beamformer  response  us 
ing  Chapter  7  structured  covariance  es 


In  Figures  30  and  31.  the  true  covariance  matrix  is  not  within  the  polyhedral  cone  formed 
by  the  fixed  basis.  There  is  a  single  interference  source  that  is  located  at  57.3  degrees  with  an 
intensity  of  60  dB  more  than  the  receiver  noise  (environment  D).  The  fixed-basis  estimator  of 
Chapter  6  shows  a  signal-to-noise  ratio  loss,  which  is  much  greater  than  the  loss  for  the  estimator 
whose  bases  are  allowed  to  vary. 

These  figures  show  that  there  can  be  an  increase  in  the  signal-to-noise  ratio  if  there  is  an 
adequate  number  of  terms  in  the  finite  sum  used  to  form  the  covariance  matrix  when  using  the 
method  of  Chapter  6;  however,  when  there  is  an  insufficient  number  of  terms,  then  this  method 
may  perform  more  poorly  than  the  unconstrained  approach.  There  is  virtually  no  change  in  the 
loss  for  the  method  of  Chapter  7  for  the  two  different  interference  environments.  This  shows  that 
the  method  of  Chapter  7  will  provide  covariance  estimates  that  are  closer  to  the  true  covariance 
matrix  for  this  case. 

9.2.3  Bias  and  Variance 

In  this  section,  the  bias  and  variance  of  the  estimators  for  the  4-element  array  are  discussed. 
It  is  known  that  all  of  the  beamforming  methods,  except  the  joint  estimation  of  the  mean  and 
covariance,  yield  unbiased  estimators  of  the  mean.  In  this  section,  simulations  are  used  to  compare 
the  bias  of  this  method  to  the  bias  of  the  other  methods.  In  this  simulation.  200  realizations 
of  the  secondaries  are  used  to  form  weight  vectors  for  the  unconstrained  and  both  constrained 
approaches;  then,  an  additional  vector  realization  of  the  same  process  is  generated.  Different  signal 
amplitudes  corresponding  to  a  desired  signal-to-noise  ratio  are  added  to  this  vector.  The  mean  is 
estimated  using  the  weight  vectors  already  formed  as  well  as  by  jointly  estimating  the  mean  with 
the  structured  covariance  matrix  using  the  estimation  methods  of  Chapters  6  and  7.  The  true  mean 
and  the  estimates  are  saved  to  a  file  where  the  sample  statistics  are  formed. 

The  covariances  that  are  used  in  these  simulations  are  within  the  polyhedral  cone  formed 
by  the  finite  basis  estimator  of  Chapter  6.  The  point  at  -20  dB  signal-to-noise  ratio  is  actually 
for  no  signal  so  that  the  bias  and  variance  of  the  estimators  with  no  input  can  be  shown.  The 
variance  for  the  known  covariance  matched  filter  is  1.  and  the  variance  of  the  estimate  that  uses 
the  unconstrained  covariance  estimate  is  2.5.  Two  interference  environments  were  simulated;  the 
first  is  environment  A.  and  the  second  is  environment  C. 

It  can  be  seen  from  Figures  32  and  33  that  the  variance  of  the  estimates  of  the  mean  using  the 
structured  methods  appears  to  be  uniformly  better  than  those  where  the  structure  of  the  covariance 
is  not  taken  into  account.  All  of  these  methods  except  the  one  that  jointly  estimates  the  mean  and 
the  covariance  are  known  to  be  unbiased.  Both  of  the  estimates  that  jointly  estimate  the  mean  and 
the  covariance  show  a  variation  in  the  bias  with  changes  in  the  signal-to-noise  ratio.  Figures  34 
and  35  are  based  on  the  interference  environment  C.  It  is  not  clear  from  these  simulations  or  others 
that  are  not  shown  if  there  is  a  trend  in  the  bias  or  variance  of  the  estimate  of  the  mean.  The  true 
mean  varies  from  1  to  20  for  the  signal-to-noise  ratio  varying  from  0  to  26  dB. 
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Figure  29.  Signal-to-noise  ratio  loss 
factor.  Chapter  7  and  unconstrained 
methods. 
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Figure  31.  Signal-to-noise  ratio  loss 
factor.  Chapter  7  and  unconstrained 
methods. 
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Figure  34-  Bias  of  the  mean  esti¬ 
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Figure  35.  Variance  of  the  mean  esti 
mates. 


9.3  Adaptive  Detection 


There  have  been  many  different  detectors  discussed  in  this  report.  These  detectors,  the 
probability  of  false  alarm,  and  the  probability  of  detection  are  shown  in  Table  5.  Of  these  detectors 


TABLE  5 


Detector  Performance  Comparisons 


Detector 

Abbreviation 

Know 

Constrained? 

Generalized  LR 

MF 

R 

- 

Generalized  LR 

GLR 

- 

N 

Adaptive  MF 

AMF 

- 

N 

Generalized  LR 

SIGI 

R* 

Y 

Generalized  LR 

CGLR 

Y 

Adaptive  MF 

CAMF 

Y 

*  Known  to  a  scale  factor 

the  first  four  have  analytic  expressions  describing  the  performance.  The  last  two  must  be  simulated. 
The  last  four  of  these  detectors  have  been  proposed  or  analyzed  as  a  portion  of  this  report.  The 
false  alarm  and  detection  probabilities  will  be  estimated  by  the  use  of  Monte  Carlo  simulations. 

Since  simulations  are  used  to  illustrate  the  performance  of  the  detection  methods,  the  proba¬ 
bilities  of  false  alarm  and  detection  will  be  restricted  to  what  would  be  considered  fairly  high  levels 
so  that  the  estimates  of  these  probabilities  are  reliable.  The  necessity  of  keeping  the  probabilities  at 
high  levels  can  be  shown  by  examining  the  variance  of  the  estimates  of  the  probabilities.  Defining 
the  indicator  function  as 

{1  x  >  Q 

(223) 

0  x  <  a 

then  for  a  single  sample  and  for  the  test  statistic  density  function  /(x), 

E{/(x)}  =  f(x)dx  =  PFA  ,  (224) 
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and 


E{(/(-r))2}  =  r  f(x)  dx  =  PFA 

Ja 


(225) 


then  the  variance  is 

a2  =  PFA  -  PFA 2  .  (226) 

For  J  realizations,  thresholding  the  test  statistic  and  averaging  the  number  of  occurrences  where 
the  test  statistic  is  greater  than  the  threshold  will  give  an  unbiased  estimate  of  the  probability  of 
detection  or  probability  of  false  alarm.  The  variance  of  the  estimate  for  the  average  of  J  realizations 
is  a2  =  (1/J)(PF.4  -  PFA2).  From  this  it  can  be  seen  that  unless  a  large  number  of  samples  is 
used  to  form  an  estimate  of  the  false  alarm  or  detection  probability,  the  variance  will  be  high.  If 
Tchebvcheff's  inequality  is  used  to  indicate  the  probability  that  the  estimate  deviates  from  the  true 
probability,  the  result  for  200  realizations  is 


P(|x  —  PFA\  >  a  *  PFA)  < 


i2  PFA2  200a2 


PFA 


-  1 


(227) 


Using  this  equation,  for  a  probability  of  detection  or  false  alarm  of  0.1.  the  probability  that  any 
realization  varies  by  more  than  50  percent  from  this  probability  is  0.18. 

The  estimates  of  the  detection  probabilities  will  not  be  entirely  independent  since  for  the 
4-element  array  only  200  covariance  estimates  will  be  performed.  These  200  estimates  are  used 
to  calculate  the  test  statistic  for  16  independent  realizations  of  the  primary  vector  at  16  different 
signal-to-noise  ratios  for  each  covariance  estimate.  For  the  full  constrained  generalized  likelihood 
ratio  test,  200  realizations  of  the  secondary  data  will  be  used  along  with  the  16  independent 
realizations  of  the  primary  vector  for  each  of  the  realizations  of  the  secondary  data.  This  yields 
3200  realizations  of  the  test  statistics  for  each  of  the  detection  methods.  For  the  8-element  array 
100  realizations  of  the  secondary  data  are  used  rather  than  200.  One  hundred  iterations  of  the 
iterative  methods  are  used  to  form  the  estimates.  For  the  joint  estimates  of  the  mean  and  the 
covariance,  the  estimate  of  the  covariance  matrix  under  Ho  was  formed;  then  this  estimate  was 
used  to  initialize  the  algorithms  for  H\. 

First,  the  two  detectors  for  which  a  closed-form  expression  for  false  alarm  probabilities  have 
not  been  found  are  simulated.  The  estimated  false  alarm  probabilities  versus  the  threshold  are 
shown  in  Figures  36  and  37  for  the  constrained  adaptive  matched-filter  detector  and  for  environ¬ 
ments  A  through  C  for  the  4-element  array.  Figures  38  and  39  illustrate  the  false  alarm  probabilities 
for  the  constrained  generalized  likelihood  ratio  detectors.  The  test  statistic  that  is  used  here  for 
the  likelihood  ratio  tests  is  the  log  of  the  K  +  1  root  of  the  ratio  of  the  likelihoods  found  by  the 
methods  of  Chapters  6  and  7. 
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The  probabilities  of  false  alarm  versus  threshold  for  the  8-element  array  are  shown  in  Fig¬ 
ures  40,  41,  42,  and  43.  It  should  be  noted  that  the  probability  of  false  alarm  is  very  sensitive  to 
changes  in  the  threshold  or  the  density  of  the  test  statistic. 

Before  discussing  the  simulation  results,  the  performance  results  for  the  detectors  with  closed- 
form  expressions  for  the  detection  probabilities  are  shown.  In  Figures  44  and  45  the  probability 
of  detection  for  the  first  four  detectors  of  Table  5  is  showm.  It  can  be  seen  that  knowing  the  noise 
statistics  adds  appreciably  to  the  detector  performance  when  the  detectors  are  compared  to  the 
approaches  that  are  not  given  any  knowledge  of  the  noise  statistics. 

The  probabilities  of  detection  for  the  detectors  introduced  in  this  report  will  now  be  compared 
to  the  other  methods.  In  order  to  plot  the  probability  of  detection,  a  desired  false  alarm  probability 
was  chosen.  The  probabilities  of  false  alarm  for  the  three  realizations  of  the  simulated  methods  are 
averaged,  and  a  threshold  is  selected  for  each  detector  and  each  estimation  method.  The  desired 
probability  of  false  alarm  was  0.1;  the  probabilities  of  false  alarm  for  the  realizations  varied  from 
0.07  to  0.135  for  the  4-element  array  and  0.06  to  0.17  for  the  8-element  array.  The  probability  of 
detection  for  the  methods  where  an  analytic  expression  exists  was  calculated  for  the  desired  false 
alarm  probability.  The  probabilities  of  detection  are  shown  for  the  three  interference  environments 
discussed  earlier  for  the  two  simulated  arrays  where  the  true  covariance  is  within  the  polyhedral 
cone. 

Figures  46  through  49  show  the  simulation  results  for  the  4-element  array,  and  Figures  50 
through  53  show  the  results  for  the  8-element  array.  It  can  be  seen  that  the  probabilities  of 
detection  at  a  particular  probability  of  false  alarm  show  a  dramatic  improvement  compared  to  the 
unconstrained  adaptive  detection  methods.  The  probability  of  detect’on  for  the  detectors  using 
the  constrained  covariance  is  nearly  that  of  the  known  covariance  matched-filter  detector.  The 
additional  signal-to-noise  ratio  that  would  be  required  to  achieve  the  same  probability  of  detection 
as  the  MF  detector  is  less  than  5  dB.  Due  to  the  statistical  variation,  there  are  some  realizations 
that  estimate  a  probability  of  detection  that  would  be  higher  than  that  of  the  knowm  covariance 
test. 

In  the  detector  simulations,  the  method  of  Chapter  6  appears  to  provide  detectors  that  have  a 
higher  probability  of  detection  compared  to  those  that  use  the  estimation  method  of  Chapter  7.  This 
is  due  to  the  location  of  the  interference  sources  with  respect  to  the  directions  that  parameterize 
the  bases.  Since  there  is  a  relatively  small  number  of  bases,  (2 L),  the  volume  of  the  constraint 
space  is  smaller  for  this  method,  and  the  truth  is  interior  to  this  smaller  volume.  The  estimator 
of  Chapter  6  effectively  has  more  knowledge  of  the  interference  environment  than  the  estimator  of 
Chapter  7.  When  the  true  covariance  matrix  is  structured,  but  not  in  the  polyhedral  cone  formed  by 
the  smaller  number  of  bases,  then  detection  performance  can  suffer.  Figure  54  illustrates  this  for  the 
interference  of  environment  D.  The  thresholds  used  to  generate  these  curves  were  those  determined 
by  the  other  simulations.  The  curves  are  labeled  CAMF6  and  CAMF7  for  the  constrained  adaptive 
matched-filter  detector  and  are  labeled  CLR6  and  CLR7  for  the  constrained  likelihood  ratio  tests. 
The  probability  of  detection  does  not  vary  appreciably  when  the  variable  bases  estimator  is  used 
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Figure  38.  Probability  of  false  alarm 
vs  threshold.  Chapter  6  likelihood  ratio. 

Figure  39.  Probability  of  false  alarm 
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Figure  41.  Probability  of  false  alarm 
vs  threshold.  Chapter  7  constrained 
AMF. 
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Figure  42 ■  Probability  of  false  alarm 
vs  threshold.  Chapter  6  likelihood  ratio. 


Figure  43  Probability  of  false  alarm 
vs  threshold.  Chapter  7  likelihood  ratio. 


114 


Figure  44  Probability  of  detection  us¬ 
ing  analytical  expressions. 
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Figure  45.  Probability  of  detection  us¬ 
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Figure  46  Probability  of  detection 
from  simulations  Chapter  fi  con¬ 
strained  .4  MF 


Figure  4~i  Probability  of  detection 
from  simulations.  Chapter  7  con¬ 
strained  AMF 


to  generate  the  test  statistics;  although,  there  is  10  dB  or  more  additional  loss  for  the  fixed-bases 
methods. 

9.4  Conclusion 

In  this  chapter,  the  adaptive  signal  processing  methods  proposed  in  this  report  have  been 
compared.  The  simulations  show  that  there  can  be  a  dramatic  improvement  in  performance  by 
making  use  of  the  complete  data  model. 
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Figure  49.  Probability  of  detection 
from  simulations.  Chapter  7  likelihood 
ratio. 
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Figure  51.  Probability  of  detection 
from  simulations.  Chapter  7  con¬ 
strained  AMF. 
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Figure  52.  Probability  of  detection 
from  simulations.  Chapter  6  likelihood 
ratio. 
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Figure  53.  Probability  of  detection 
from  simulations.  Chapter  7  likelihood 
ratio. 


Figure  54.  Probability  of  detection 
from  simulations ;  Chapter  6  con¬ 
strained  AMF 


10.  CONCLUSIONS 


The  subject  of  this  report  has  been  the  use  of  sensor  arrays  to  perform  beamforming  and 
detection.  Maximum-likelihood  techniques  have  been  used  to  develop  approaches  which  can  be 
implemented  to  accomplish  these  tasks.  The  results  of  this  report  and  possible  future  research 
areas  are  summarized  below. 

The  radar  detection  and  the  communications  applications  of  this  report  were  discussed  in 
Chapter  2,  and  the  resulting  data  model  was  introduced.  The  radar  detection  model  is  the  well- 
known  non-fluctuating  target  model  for  which  target  presence  is  indicated  by  a  non-zero  structured 
mean.  The  structured  mean  is  modeled  as  an  unknown  scalar  times  an  array  steering  or  direction 
vector.  In  this  model,  the  covariance  matrix  is  composed  of  two  terms.  The  first  is  a  diagonal 
matrix  due  to  noise  sources  internal  to  the  receivers;  the  second  is  generated  by  an  integral  ex¬ 
pression  consisting  of  a  positive  spatial-temporal  power-spectral  density  and  the  outer  products  of 
the  array  direction  vectors.  The  communications  model  is  similar  to  the  detection  model,  with  the 
transmitted  information  content  contained  in  the  structured  mean. 

There  are  several  modifications  of  the  data  model  that  would  lead  to  areas  of  future  research. 
One  model  could  remove  the  assumption  that  the  direction  of  a  target  and  the  data  vector  containing 
a  possible  target  return  are  known.  Multiple  targets  could  then  be  allowed.  As  shown  in  Chapter  6, 
the  likelihood  can  be  unbounded  for  some  conditions;  increasing  the  number  of  sample  vectors  that 
can  contain  a  mean  will  make  this  problem  worse.  A  straightforward  maximization  of  the  likelihood 
may  not  be  feasible,  and  methods  of  preventing  this  from  occurring  would  need  to  be  developed. 

An  additional  model  which  could  be  the  subject  of  further  research  would  be  to  remove  the 
assumption  of  independence  for  the  data  samples.  A  detector  which  could  detect  a  moving  target 
in  clutter  based  on  the  returns  from  many  samples  could  then  be  developed. 

Chapter  3  is  used  to  review  some  of  the  current  methods  of  beamforming  and  detection.  It 
was  pointed  out  that  optimal  beamforming  and  detection  require  knowledge  of  the  noise  statistics. 
As  these  statistics  are  seldom  known,  adaptive  approaches  are  utilized  to  estimate  the  statistics.  If 
the  amount  of  data  is  limited,  then  there  can  be  an  appreciable  loss  in  performance  for  the  adaptive 
approaches  motivating  the  use  of  a  more  detailed  data  model. 

In  Chapter  4,  a  detector  is  derived  and  analyzed.  The  derivation  uses  the  structured  mean  of 
the  data  model  but  does  not  use  any  knowledge  of  the  interference  environment  or  the  array  when 
estimating  the  noise  covariance  matrix.  This  detector  provides  a  desirable  simplification  of  the  full 
generalized  likelihood  ratio  test  and  provides  similar  performance  for  signals  which  are  in  alignment 
with  the  assumed  signal  direction.  An  original  contribution  discussed  in  this  report  is  the  analysis 
of  detection  performance  when  the  signal  is  not  aligned  with  the  assumed  signal  direction. 

The  set  of  noise  covariance  matrices  which  is  possible  for  a  particular  array  under  the  model 
discussed  in  Chapter  2  is  the  subject  of  Chapter  5.  The  set  of  possible  covariance  matrices  is 
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characterized  as  a  convex  cone  that  is  a  subset  of  an  N 2  dimensional  real  vector  space.  Struc¬ 
tured  covariance  matrices  are  defined  to  be  a  convex  cone  in  a  lower  dimensional  space.  Some 
of  the  representations  for  covariance  matrices  that  are  structured  are  introduced  in  this  chapter. 
These  representations  are  introduced  with  the  intent  that  they  be  used  to  enforce  the  structure  on 
estimated  covariance  matrices  that  are  members  of  the  constraint  space. 

One  of  the  representations  introduced  in  Chapter  5  approximates  the  integral  expression  for 
the  constraint  space  by  a  finite  sum.  The  result  consists  of  a  sum  of  outer  products  of  the  array 
direction  vectors  weighted  by  samples  of  the  spatial-temporal  power-spectral-density  function.  A 
joint  maximum-likelihood  estimator  of  the  structured  mean  and  covariance  matrix  are  developed  in 
Chapter  6.  The  resulting  iterative  procedure  is  based  on  the  Expectation-Maximization  algorithm 
and  results  in  a  sequence  of  estimates  for  which  the  log- likelihood  is  non-decreasing.  Stable  points 
of  the  resulting  estimator  satisfy  the  necessary  conditions  for  a  maximizer  of  the  likelihood  function. 
Certain  technical  aspects  of  the  resulting  estimation  algorithm  are  discussed  in  this  chapter  as  well. 

Properties  of  the  likelihood  surface  and  a  maximum-likelihood  estimate  of  the  weight  of  single 
spectral  point  are  discussed  as  an  introduction  to  Chapter  7.  Based  on  this  discussion  and  one  of 
the  representations  of  Chapter  5,  an  algorithm  for  estimating  constrained  means  and  covariance 
matrices  is  introduced.  This  method  does  not  require  the  large  number  of  spectral  weights  but 
provides  resolution  and  performance  advantages  compared  to  the  estimator  of  Chapter  6.  The 
likelihood  for  the  sequence  of  estimates  produced  by  this  estimator  is  non-decreasing,  and  stable 
points  of  the  resulting  algorithm  satisfy  the  necessary  conditions  for  a  maximizer  of  the  likelihood. 

Applying  the  structured  estimators  to  beamforming  and  detection  is  the  subject  of  Chapter  8. 
Two  methods  of  utilizing  the  estimates  produced  in  Chapters  6  and  7  to  perform  both  beamforming 
and  detection  are  discussed.  Several  properties  of  these  methods  are  shown,  including  bias  and  some 
aspects  of  CFAR  behavior. 

The  methods  of  beamforming  and  detection  discussed  in  this  report  are  compared  in  Chap¬ 
ter  9.  It  is  shown  that  there  can  be  a  dramatic  increase  in  signal-to-noise  ratio  for  the  adaptive 
beamformers  based  on  the  structured  covariance  estimates  compared  to  beamformers  which  are  not 
based  on  this  knowledge.  This  is  seen  in  the  signal-to-noise  ratio  loss  factor  and  in  the  beamformer 
responses.  For  the  detectors  which  use  knowledge  of  the  array  geometry  and  the  noise  environment, 
there  is  a  similar  increase  in  performance.  The  false  alarm  probability  estimated  by  the  simulations 
has  a  statistical  variation  that  would  be  within  what  would  be  expected  if  the  detectors  had  the 
CFAR  property.  The  probability  of  detection  shows  the  dramatic  increase  in  performance  that 
was  seen  in  the  constrained  adaptive  beamformers.  The  loss  in  signal-to-noise  ratio  for  the  same 
probability  of  detection  is  much  less  than  the  methods  which  do  not  utilize  the  knowledge  of  the 
covariance  structure. 
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APPENDIX  A 

a2 1  CONSTRAINED  DETECTOR 


A  convenient  test  that  may  be  used  as  a  reference  occurs  when  the  covariance  is  constrained 
to  be  of  the  form  o2I.  When  the  covariance  is  known  with  the  exception  of  a  scale  factor,  the 
results  of  this  appendix  may  be  used  by  transforming  the  data  to  the  cr2I  case  through  whitening. 

The  notation  used  here  is  identical  to  that  used  in  Chapter  4.  A  single  primary  vector  z  is 
zero  mean  on  Ho  and  has  a  mean  of  the  form  bd  under  H\ .  In  addition,  K  mutually  independent 
zero  mean  secondary  vectors,  the  columns  of  Z,  are  available. 

For  a  single  vector  sample  with  mean  bd ,  the  complex  Gaussian  density  function  is 

/( z)  =  7r-va-2've^(2-bd)t(z-6d)  (A.l) 


The  generalized  likelihood  ratio  test  over  the  joint  density  function  of  (z.Z)  is 


ma*  fz.Z'Hi  (z-  b,  a2\H\) 

\  = 


max/z.Z|//0(z-Z;aoltfo) 


^  a 

Ho 


(A.2) 


Substituting  in  the  density  function  and  canceling  common  terms  yields 


max  a 

A=±l 


1  e  1 


max(7c 


2,V( A'+l)  ^!(z’z~£.=  i 


(A- 3) 


e  o 


The  maximum  likelihood  estimate  for  o\  is  found  by  maximizing  the  numerator  with  respect 
to  of  with  the  result  that 


°\  =  iv(^TT)^Z_  bd)\z-bd)  +  J2*l*i}  >and 


(A.4) 


likewise 


On 


i 


K 


N(K  +  1) 


z'z  + 


,  zIzl 


i=l 


(A.5) 
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Substituting  these  estimates  into  Equation  (A. 3)  yields 


A  = 


,v(a-+T)[(z  "  fcd)t(z  —  6d)  +  N' 


K-<-l 


N{K+\)  1 


■N(K  + 1) 
e-X(K+l) 


(A-6) 


Canceling  common  terms  and  taking  the  N(K  +  1)  root  (a  monotonic  functional  of  the  likelihood) 
yields 


A  = 


max[(z  -  6d)t(z  -  6d)  +  z!zii  1  Hi 


[zfz  +  Etil  Az' 


Ho 


a 


(A. 7) 


a  has  been  redefined  to  include  the  effects  of  the  N(K  +  1)  root. 

Maximizing  the  numerator  with  respect  to  b  is  equivalent  to  minimizing 

.  *  Id^zi2  +  d*z  o 

(z  -  6d)f(z  -  6d)  =  zfz  -  +  dtd|-^  -  b\2 

which  is  clearly  minimized  by 


Substituting  Equation  (A. 9)  into  Equation  (A. 7)  yields 


(A- 8) 


(A.9) 


A  = 


+  £,=i  z!z.  ^ 


zTz  -  +  E«  =  l  z,zi 


^  a 

Ho 


(A.10) 


Or  rewriting  yields 


|dtz!2 

dtd 


Hi 

Ho 


0-1 

a 


\z^z 


K 

1=1 


ZIZ«1 


(All) 
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Thus,  we  are  comparing  the  power  of  the  estimate  of  6  with  the  sum  of  the  power  of  the  other 
vector  components,  which  we  will  see  is  a  simple  scalar  CFAR  detector. 

The  vector  notation  for  the  decision  rule  can  now  be  rewritten  where  the  test  will  be  more 
readily  recognized  as  a  scalar  CFAR  test,  and  the  performance  may  be  found. 

Perform  a  unitary  transform  on  the  data  of  the  form 

y  =  Ufz  U  =  [d:Uj]  ,  (A.12) 

and  Y  =  tPZ.  The  normalizing  assumption  will  be  made  that  dM  =  1.  The  test  can  be 
rewritten  in  this  notation  as 


^UU^i2 


!dfUy! 


H , 

£ 

H0 


a  —  1 


Q 


K 


iyfy  +  I>,ty.! 


i=l 


(A. 13) 


cFU  is  the  elementary  vector  [1,0,  ...,0j,  so  the  test  can  be  written  in  terms  of  the  vector  elements 
as 


!yii2  ^ 


Hy  _  _  1  x  K  .V 

'  ‘■Ew’  +  LDn.il 


Q 


H0  1=1  »=lj=l 


or 


"i  (a-DE^  +  EtlVl 

H0  1=2  i=ll=l 


(A. 14) 


(A. 15) 


This  form  of  the  detector  is  now  easily  recognized  as  a  scalar  CFAR  detector.  The  data  can 
be  normalized  so  that  the  variance  of  the  y/s  is  unity.  Then  the  y/s  are  then  distributed  N( 0, 1), 
and  y\  is  distributed  N( 0, 1)  or  N(^,  1). 

Under  Hq,  |yi|2  is  distributed  chi-square  with  2  degrees  of  freedom,  Y.^2  i J/j I2  is  distributed 
chi-square  with  2 (N  -  1)  degrees  of  freedom,  and  H)V=i  ! Yj  1 2  distributed  chi-square  with  2 KN 

degrees  of  freedom.  The  right-hand  side  of  Equation  (A. 15)  has  2[(A’  +  l)Ar  -  lj  degrees  of  freedom, 
and  we  may  use  these  distributions  to  evaluate  the  performance.  The  probability  of  false  alarm  is 
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noo  t1-  1  i 

-yj-  e~T  e~w  dwdT  =  -j 
a-l)  L!  Q 


(A. 16) 


where  L  =  (if  +  1  )N  —  1 

The  probability  of  detection  is  found  through  comparison  with  Kelly  [32]: 


PDgigI  —  1  r 


4  if  M(a-DfcC*(-),  Gk(y)  =  e-ykfy-  ,  (A.17) 

aLt i\k  )  „=on! 


where  a  is  the  signal-to- noise  ratio  |fe|2/ cr2. 

If  the  covariance  estimate  based  on  this  model  is  used  to  perform  adaptive  beamforming,  then 
there  will  be  no  loss  in  signal- to- noise  ratio.  Looking  at  the  loss  factor  density, 


IdtRMl2 

dtR-MdtR-iRR-M 


(A. 18) 


and  substituting  cr2 1  for  R,  the  result  is 


[dfd|2 

P  ~  dtR  iddtRd 


(A.19) 


This  is  independent  of  the  estimate  of  a2  and  is  unity  provided  that  the  model  is  matched. 
When  the  unconstrained  maximum-likelihood  estimate  of  the  covariance  matrix  was  used  to  perform 
beamforming  the  loss  factor  is  a  random  variable  [10].  For  this  case,  the  loss  factor  is  not  a  random 
variable. 
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GLOSSARY 


N 

G 

K 

P 

M 

L 


0, 0 

S(lj,6,0) 

7 

A 

l 

Z.  (*n) 

z 

B 

b 

R 

R 

s 

d,p,q 

D 

£ 

L 

n, 

a 


Number  of  sensors,  data  vector  length 
Number  of  non-zero  mean  data  vectors 
Number  of  zero  mean  data  vectors 
Number  of  data  vectors  (K  +  G) 

Number  of  terms  in  finite  sum  approximation 
Dimension  of  covariance  matrix  vector  space 
Sample  index 
Carrier  frequency 

Physical  angles  relative  to  array  reference 

Spatial-temporal  power-spectral  density 

Threshold 

Likelihood  ratio 

Log-likelihood 

Sampled  data  vector  (at  time  n) 

Data  matrix  containing  column  vectors  z\  . .  .zp 
Matrix  which  is  E(Z) 

Scalar  portion  of  a  structured  mean,  E(z)  =  bd 

Covariance  matrix 

Estimated  covariance  matrix 

Array  coordinate  vector 

Array  steering  vector 

Matrix  whose  columns  are  array  steering  vectors 

Diagonal  matrix  of  discrete  spectral  weights 

Transformation,  usually  a  unitary  matrix 

The  space  of  possible  covariance  matrices  for 
a  particular  array  geometry 

Signal-to-noise  ratio  available  using  optimal  beamformer 
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GLOSSARY 

(Continued) 

Ratio  of  signal-to-noise  ratio  for  a  linear  beamformer  to 
optimal  signal-to-noise  ratio 

Real  value  of  argument 

Imaginary  value  of  argument 
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adaptive  beamforming  and  detection.  A  constant  false  alarm  rate  detector,  based  on  unconstrained  maximum-likelihood 
techniques,  is  derived  and  analyzed.  The  spare  of  noise  covariance  matrices  possible  from  a  particular  array  is  character¬ 
ized.  vielding  representations  for  the  spare  and  members  of  the  space  in  terms  of  finite  numbers  of  spectral  points.  These 
representations  are  used  to  derive  constrained  maximum-likelihood  estimators  that  jointly  estimate  the  parameters  of  the 
density  function.  Two  approaches  that  use  the  constrained  covariance  estimates  to  perform  beamfonning  are  described  and 
compared.  The  loss  in  signal-lo-noise  ratio  and  the  v  ariance  of  the  estimators  are  shown  to  be  less  for  these  approaches  than 
for  those  that  do  not  use  the  covariance  model.  Detection  methods  based  on  the  generalized  likelihood  ratio  test  and  a 
constant  false  alarm  rate  matched-filter  detector  are  analvzed.  and  simulation  results  are  presented. 
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